text:
"""
   FOR PRE PROCESSING ONLY (generates files in test/data)
"""
function _genpolmeshes(; numb::Int64 = 2,
                         nbit::Union{Vector{Int64}, Int64} = [100,],
                         np::Union{Vector{Int64}, Int64} = [100, 1_000, 10_000],
                         drawvoronoi::Bool = false,
                         dim::Int64 = 2)
    support = if dim == 1 && numb == 2
        "disk"
    elseif dim == 2 && numb == 2
        "square"
    elseif dim == 2 && numb == 4
        "flower"
    elseif dim == 2 && numb == 9
        "L"
    elseif dim == 3 && numb == 1
        "bar"
    elseif dim == 3 && numb == 3
        "cube"
    elseif dim == 3 && numb == 6
        "ico"
    end
    # gen regular polygonal meshes using Lloyd
    for itmax ∈ nbit
        strit = if itmax == 1
            "_coarse"
        elseif itmax < 1_000
            "_lesscoarse"
        else
             ""
        end
        for npl ∈ np
            # WARNING function not available (only for local use)
            c = Main.ConvexTools.lloyd_geogram(npl, numb = numb, nbit = itmax, dim = dim,
                          drawvoronoi = drawvoronoi)
            # save
            mesh_filename = "$(@__DIR__)/../test/data/$(support)polmesh$(strit)_$(npl).jld2"
            println(mesh_filename)
            JLD2.save(mesh_filename, "ps", c.ps, "pv", c.pv, "elem", c.elem,
                      "bbelem", c.bbelem, "pb", c.pb, "tb", c.boxelem,
                      "Iv", c.Iv, "centroids", c.centroids, "measures", c.measures)
        end
    end
    nothing
end
function _genallpolmeshes()
    nbit = [1, 100, 1_000]
    np = [100, 1_000, 10_000]
    dimnum = [(2, 2), (2, 4), (2, 9), (3, 1), (3, 3), (3, 6)]
    for it ∈ nbit
        for nn ∈ np
            for v ∈ dimnum
                dim, numb = v
                _genpolmeshes(numb = numb, dim = dim, nbit = it,
                              drawvoronoi = false, np = nn)
            end
        end
    end
    nothing
end
#-------------------------------------------------------------------------------
"""
   FOR PRE PROCESSING ONLY (generates files in test/data)

CALL: _genmeshes(numb = 2, np = [100, 1_000, 10_000])
"""
function _genmeshes(; numb::Int64 = 2,
                      np::Union{Vector{Int64}, Int64} = [100, 1_000, 10_000],
                      drawvoronoi::Bool = false)
    support = if numb == 2
        "square"
    elseif numb == 4
        "flower"
    elseif numb == 9
        "L"
    end
    # gen regular meshes calling triangle
    for npl ∈ np
        # WARNING function not available (only for local use)
        _, ptri, tri = Main.MeshTools.triangle(numb, nbtri = 2 * npl)
        # keep only interior points for voronoi cells
        Ik = setdiff(1:size(ptri, 1), unique(btri(tri)[:]))
        # WARNING function not available (only for local use)
        p, pv, cellsb, cellsbt, t, pc, pb, tb = Main.ConvexTools.lloyd_geogram(npl, numb = numb, nbit = 1, p = copy(ptri[Ik, :]), drawvoronoi = drawvoronoi)
        if numb == 2 # prefere a unit uncentered square
            ptri .+= 1
            ptri ./= 2
        end
        # save
        mesh_filename = "$(@__DIR__)/../test/data/$(support)mesh_$(npl).jld2"
        JLD2.save(mesh_filename, "pv", pv, "cellsb", cellsb,
                  "cellsbt", cellsbt, "pb", pb, "tb", tb, "t", t,
                  "ptri", ptri, "tri", tri)
        println(mesh_filename)
    end
    nothing
end
#-------------------------------------------------------------------------------
"""
   FOR POST PROCESSING ONLY
"""
function _genfigures()
    for file ∈ ["Lpolmesh", "Lpolmesh_lesscoarse", "Lpolmesh_coarse",
                "squarepolmesh", "squarepolmesh_lesscoarse",
                "squarepolmesh_coarse"]
        for np ∈ [100, 1000, 10000]
            mesh_filename = "$(@__DIR__)/../test/data/$(file)_$(np).jld2"
            vem(file, np, resolution = 1_000)
            figfile = replace(replace(mesh_filename, "jld2" => "png"),
                                                     "data" => "figures")
            println("save figure to $(figfile)")
            Makie.save(figfile, Makie.current_figure())
        end
    end
    nothing
end
#-------------------------------------------------------------------------------
"""
   FOR POST PROCESSING ONLY

CALL: _gensolutions(); _generrors(); _plotorder("square")
"""
function _gensolutions()
    sols = Dict()
      # P1 solutions
      for file ∈ ["Lmesh", "squaremesh"]
        for np ∈ [100, 1000, 10000]
            mesh_filename = "$(@__DIR__)/../test/data/$(file)_$(np).jld2"
            u, ptri, tri, Ib = P1(file, np, resolution = 0)
            sols[mesh_filename] = (u, ptri, tri, Ib)
            mesh_filename = replace(mesh_filename, "mesh" => "polmeshP1")
            u, pv, t, Ib = vem(file, np, resolution = 0)
            sols[mesh_filename] = (u, pv, tri, Ib)
        end
    end
    # Virtual method solutions
    for file ∈ ["Lpolmesh", "Lpolmesh_lesscoarse", "Lpolmesh_coarse",
                "squarepolmesh", "squarepolmesh_lesscoarse",
                "squarepolmesh_coarse"]
        for np ∈ [100, 1000, 10000]
            mesh_filename = "$(@__DIR__)/../test/data/$(file)_$(np).jld2"
            u, pv, t, Ib = vem(file, np, resolution = 0)
            sols[mesh_filename] = (u, pv, t, Ib)
            println("solution of $(mesh_filename)")
        end
    end
    @save "/tmp/sols.jld2" sols
    nothing
end
#-------------------------------------------------------------------------------
"""
   FOR POST PROCESSING ONLY
"""
function _generrors()
    @load "/tmp/sols.jld2"
    fsolsquare = "$(@__DIR__)/../test/data/squaresolution_100000.jld2"
    fsolL = "$(@__DIR__)/../test/data/Lsolution_100000.jld2"
    dsolsquare, dsolL =  FileIO.load(fsolsquare), FileIO.load(fsolL)
    errors = Dict()
    for ff ∈ keys(sols)
        println(ff)
        dsol = if occursin("square", ff)
            dsolsquare
        else
            dsolL
        end
        # compute error in ||.||∞ for inside points
        Inb = setdiff(1:size(sols[ff][2], 1), sols[ff][4])

        Itri, Iclose, baryc = Main.barycentric_coordinates(dsol["ptri"],
                                                           dsol["tri"],
                                                           sols[ff][2])
        Inb = findall(Itri .> 0)
        println("percentage of interior points ", length(Inb) / size(Itri, 1) *
                                                  100)
        ui = sum(dsol["u"][dsol["tri"][Itri[Inb], k]] .* baryc[Inb, k] for k = 1:3)
        err = maximum(abs.(ui .- sols[ff][1][Inb]))
        err = norm(ui .- sols[ff][1][Inb]) / sqrt(length(Inb))
        nbdof = length(sols[ff][1])
        errors[ff] = (err, nbdof)
    end
    @save "/tmp/errors.jld2" errors
    nothing
end
function _ploterror(file::String = "squarepolmesh_coarse", np::Int64 = 1_00)
    @load "/tmp/sols.jld2"
    mesh_filename = "$(@__DIR__)/../test/data/$(file)_$(np).jld2"
    ff = mesh_filename
    mesh_filename = replace(mesh_filename, "polmeshP1" => "mesh")
    fsolsquare = "$(@__DIR__)/../test/data/squaresolution_100000.jld2"
    fsolL = "$(@__DIR__)/../test/data/Lsolution_100000.jld2"
    dsolsquare, dsolL =  FileIO.load(fsolsquare), FileIO.load(fsolL)
    dsol = if occursin("square", ff)
        dsolsquare
    else
        dsolL
    end
    # compute error in ||.||∞ for inside points
    Inb = setdiff(1:size(sols[ff][2], 1), sols[ff][4])
    Itri, Iclose, baryc = Main.barycentric_coordinates(dsol["ptri"],
                                                       dsol["tri"],
                                                       sols[ff][2])
    Inb = findall(Itri .> 0)
    println("percentage of interior points ", length(Inb) / size(Itri, 1) *
                                              100)
    ui = sum(dsol["u"][dsol["tri"][Itri[Inb], k]] .* baryc[Inb, k] for k = 1:3)
    err = zeros(size(Itri, 1))
    err[Inb] = abs.(ui .- sols[ff][1][Inb])
    if !occursin("pol", file)
        GLMakie.destroy!(GLMakie.global_gl_screen())
        resolution = 400
        JLD2.@load(mesh_filename, ptri, tri)
        f, L = Main.figure(1, fig3D = false, resolution = (resolution, resolution))
        Main.plot_t(ptri, tri, representation = "meshsurface",
                    scalars = err, colormap = "hsv", colorw = (0., 0, 0),
                    scene = L[1])
        Colorbar(f[1, 2], limits = (minimum(err), maximum(err)),
                 colormap = :hsv,
                 size = 10,  height = Relative(3/4))
    else
        JLD2.@load(mesh_filename, pv, cellsb)
        plotsolution(err, pv, cellsb, mesh_filename = mesh_filename,
                     resolution = 400)
    end
    nothing
end
#-------------------------------------------------------------------------------
"""
   FOR POST PROCESSING ONLY
"""
function _plotorder(geom::String = "square")
    @load "/tmp/errors.jld2"
    GLMakie.destroy!(GLMakie.global_gl_screen())
    fig = GLMakie.Figure(resolution = (1200, 400))
    ax = L = GLMakie.Axis(fig[1, 1], xlabel = "log nb dof", ylabel = "log L² error ")
     plotsl, plotsm = Any[], Any[]
    markers = [:circle, :rect, :utriangle, :dtriangle, :diamond]
    allfiles = if occursin("square", geom)
        ["squarepolmesh", "squarepolmesh_lesscoarse", "squarepolmesh_coarse",
          "squaremesh",  "squarepolmeshP1"]
    else
        ["Lpolmesh", "Lpolmesh_lesscoarse", "Lpolmesh_coarse", "Lmesh",
         "LpolmeshP1"]
    end
    cm = 1
    colors = [RGB(0.8, 0.2, 0.1), RGB(0.8, 0.8, 0.1), RGB(0.1, 0.8, 0.2),
              RGB(0.1, 0.1, 0.8), RGB(0.1, 0.7, 0.8)]
    for file ∈ allfiles
        x, y = Float64[], Float64[]
        for np ∈ [100, 1000, 10000]
            mesh_filename = "$(@__DIR__)/../test/data/$(file)_$(np).jld2"
            err, nbdof = errors[mesh_filename]
            push!(x, log(nbdof))
            push!(y, log(err))

        end
        cr = colors[cm] #rand(RGB)
        push!(plotsl, Makie.lines!(x, y, linewidth = 5, color = cr))
        push!(plotsm, Makie.scatter!(x, y, color = cr, marker = markers[cm],
              markersize = 12))
        cm += 1
    end
    ax = Axis(fig[1, 1])
    Makie.Legend(fig[2, 1], [[plotsl[k], plotsm[k]] for k = 1:length(allfiles)],
                 allfiles, tellwidth = false, tellheight = true,
                 orientation = :horizontal)
    display(Makie.current_figure())
    figfile = if occursin("square", geom)
        "$(@__DIR__)/../test/figures/convergence_square.png"
    else
        "$(@__DIR__)/../test/figures/convergence_L.png"
    end
    println("save figure to $(figfile)")
    Makie.save(figfile, Makie.current_figure())
    display(Makie.current_figure())
    nothing
end


text:
using NetworkEpidemics
using LightGraphs
using Plots

N = 1000
g = complete_graph(N)
g = erdos_renyi(N, 0.3)
g = random_regular_graph(N,15)
M = 10
h = complete_graph(M)

β = 0.1
D = 0.1

# SI

mp = Metapopulation(h, D, SI(β))
mpx = Metaplex(g, h, D, SI(β))

x0_mp = Array{Int,2}(undef, M,2)
x0_mp[2:M,:] .= [100 0]
x0_mp[1, :] .= [90, 10]

x0_i = fill(1, N)
x0_i[1:10] .= 2
x0_μ = [div(i-1,div(N,M))+1 for i in 1:N]

tmax = 100.0
nmax = 12000

@time ts_mp, output_mp = gillespie(mp, x0_mp, nmax=nmax)
@time ts_av_mp, u_av_mp = average(mp, x0_mp, nmax=nmax, tmax=tmax, progressbar=true)
@time ts_mf_mp, u_mf_mp  = meanfield(mp, x0_mp, tmax=tmax, saveat=ts_av_mp)

@time ts_mpx, output_mpx = gillespie(mpx, [x0_i,x0_μ], nmax=nmax)
@time ts_av_mpx, u_av_mpx = average(mpx, [x0_i,x0_μ], nmax=nmax, tmax=tmax, progressbar=true)
@time ts_mf_mpx, u_mf_mpx = meanfield(mpx, [x0_i, x0_μ], tmax=tmax, saveat=ts_av_mpx)

plot(ts_mp, output_mp[2],
    line=:steppre,
    label="",
    title="Metapopulation")
plot(ts_mpx, output_mpx[2],
    line=:steppre,
    label="",
    title="Metaplex")

plot(ts_av_mp, u_av_mp[2],
    label="",
    linestyle=:dash,
    title="Metapopulation")
plot!(ts_mf_mp, u_mf_mp[2], label="")

plot(ts_av_mpx, u_av_mpx[2],
    label="",
    linestyle=:dash,
    title="Metaplex")
plot!(ts_mf_mpx, u_mf_mpx[2], label="")

plot(ts_av_mp, sum(u_av_mp[2], dims=2),
    label="Average (Metapopulation)",
    linestyle=:dash,
    legend=:bottomright,
    title="Metapopulation vs Metaplex (SI) - RRG(k=15)"
    )
plot!(ts_mf_mp, sum(u_mf_mp[2], dims=2),
    label="Meanfield (Metapopulation)")
plot!(ts_av_mpx, sum(u_av_mpx[2], dims=2),
    label="Average (Metaplex)",
    linestyle=:dash,
    )
plot!(ts_mf_mpx, sum(u_mf_mpx[2], dims=2),
    label="Meanfield (Metaplex)")

# SIS
β = 0.1
γ = 1.0
D = 0.1

mp = Metapopulation(h, D, SIS(β,γ))
mpx = Metaplex(g, h, D, SIS(β,γ))
c_p = ContactProcess(g, SIS(β,γ))

x0_mp = Array{Int,2}(undef, M,2)
x0_mp[2:M,:] .= [100 0]
x0_mp[1, :] .= [90, 10]

x0_i = fill(1, N)
x0_i[1:10] .= 2
x0_μ = [div(i-1,div(N,M))+1 for i in 1:N]

x0_cp = x0_i .== 2

tmax = 10.0
nmax = 8000

@time ts_mp, output_mp = gillespie(mp, x0_mp, nmax=nmax)
@time ts_av_mp, u_av_mp = average(mp, x0_mp, nmax=nmax, tmax=tmax, progressbar=true)
@time ts_mf_mp, u_mf_mp  = meanfield(mp, x0_mp, tmax=tmax, saveat=ts_av_mp)

@time ts_mpx, output_mpx = gillespie(mpx, [x0_i,x0_μ], tmax=tmax, nmax=nmax)
@time ts_av_mpx, u_av_mpx = average(mpx, [x0_i,x0_μ], nmax=nmax, tmax=tmax, progressbar=true)
@time ts_mf_mpx, u_mf_mpx = meanfield(mpx, [x0_i, x0_μ], tmax=tmax, saveat=ts_av_mpx)

@time ts_cp, output_cp = gillespie(c_p, x0_cp, nmax=nmax)
@time ts_av_cp, u_av_cp = average(c_p, x0_cp, nmax=nmax, tmax=tmax, progressbar=true)
@time ts_mf_cp, u_mf_cp = meanfield(c_p, x0_cp, tmax=tmax, saveat=ts_av_cp)

plot(ts_mp, output_mp[2],
    line=:steppre,
    label="",
    title="Metapopulation")
plot(ts_mpx, output_mpx[2],
    line=:steppre,
    label="",
    title="Metaplex")
plot(ts_cp, output_cp,
    line=:steppre,
    label="",
    title="Contact Process")

plot(ts_av_mp, u_av_mp[2],
    label="",
    linestyle=:dash,
    title="Metapopulation")
plot!(ts_mf_mp, u_mf_mp[2], label="")

plot(ts_av_mpx, u_av_mpx[2],
    label="",
    linestyle=:dash,
    title="Metaplex")
plot!(ts_mf_mpx, u_mf_mpx[2], label="")

plot(ts_av_mp, sum(u_av_mp[2], dims=2)/N,
    label="Average (Metapopulation)",
    linestyle=:dash,
    legend=:bottomright,
    title="Metapopulation vs Metaplex (SIS) - RRG(k=15)",
    xlabel="Time",
    ylabel="Fraction of infected individuals"
    )
plot!(ts_mf_mp, sum(u_mf_mp[2], dims=2)/N,
    label="Meanfield (Metapopulation)")
plot!(ts_av_mpx, sum(u_av_mpx[2], dims=2)/N,
    label="Average (Metaplex)",
    linestyle=:dash,
    )
plot!(ts_mf_mpx, sum(u_mf_mpx[2], dims=2)/N,
    label="Meanfield (Metaplex)")
plot!(ts_mf_cp, sum(u_mf_cp, dims=2)/N,
    label="Meanfield (Contact Process)")

# SIR
β = 0.1
γ = 1.0
D = 0.1

mp = Metapopulation(h, D, SIR(β,γ))
mpx = Metaplex(g, h, D, SIR(β,γ))

x0_mp = Array{Int,2}(undef, M, 3)
x0_mp[2:M,:] .= [100 0 0]
x0_mp[1, :] .= [90, 10, 0]

x0_i = fill(1, N)
x0_i[1:10] .= 2
x0_μ = [div(i-1,div(N,M))+1 for i in 1:N]

tmax = 18.0
nmax = 5000

@time ts_mp, output_mp = gillespie(mp, x0_mp, nmax=nmax)
@time ts_av_mp, u_av_mp = average(mp, x0_mp, nmax=nmax, tmax=tmax, progressbar=true)
@time ts_mf_mp, u_mf_mp  = meanfield(mp, x0_mp, tmax=tmax, saveat=ts_av_mp)

@time ts_mpx, output_mpx = gillespie(mpx, [x0_i,x0_μ], tmax=tmax, nmax=nmax)
@time ts_av_mpx, u_av_mpx = average(mpx, [x0_i,x0_μ], nmax=nmax, tmax=tmax, progressbar=true)
@time ts_mf_mpx, u_mf_mpx = meanfield(mpx, [x0_i, x0_μ], tmax=tmax, saveat=ts_av_mpx)

plot(ts_mp, output_mp[2],
    line=:steppre,
    label="",
    title="Metapopulation")
plot(ts_mpx, output_mpx[2],
    line=:steppre,
    label="",
    title="Metaplex")

plot(ts_av_mp, u_av_mp[2],
    label="",
    linestyle=:dash,
    title="Metapopulation")
plot!(ts_mf_mp, u_mf_mp[2], label="")

plot(ts_av_mpx, u_av_mpx[2],
    label="",
    linestyle=:dash,
    title="Metaplex")
plot!(ts_mf_mpx, u_mf_mpx[2], label="")

plot(ts_av_mp, sum(u_av_mp[2], dims=2)/N,
    label="Average (Metapopulation)",
    linestyle=:dash,
    legend=:topright,
    title="Metapopulation vs Metaplex (SIR) - RRG(k=15)",
    xlabel="Time",
    ylabel="Fraction of infected individuals"
    )
plot!(ts_mf_mp, sum(u_mf_mp[2], dims=2)/N,
    label="Meanfield (Metapopulation)")
plot!(ts_av_mpx, sum(u_av_mpx[2], dims=2)/N,
    label="Average (Metaplex)",
    linestyle=:dash,
    )
plot!(ts_mf_mpx, sum(u_mf_mpx[2], dims=2)/N,
    label="Meanfield (Metaplex)")

plot(ts_av_mp, sum(u_av_mp[3], dims=2)/N,
    label="Average (Metapopulation)",
    linestyle=:dash,
    legend=:bottomright,
    title="Metapopulation vs Metaplex (SIR) - RRG(k=15)",
    xlabel="Time",
    ylabel="Fraction of Recovered individuals"
    )
plot!(ts_mf_mp, sum(u_mf_mp[3], dims=2)/N,
    label="Meanfield (Metapopulation)")
plot!(ts_av_mpx, sum(u_av_mpx[3], dims=2)/N,
    label="Average (Metaplex)",
    linestyle=:dash,
    )
plot!(ts_mf_mpx, sum(u_mf_mpx[3], dims=2)/N,
    label="Meanfield (Metaplex)")

ext: jl

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
using PowerSystems
using PowerSimulations
using JuMP
using Cbc

base_dir = string(dirname(dirname(pathof(PowerSystems))))
println(joinpath(base_dir,"data/data_5bus.jl"))
include(joinpath(base_dir,"data/data_5bus.jl"))

sys5b = PowerSystem(nodes5, generators5, loads5_DA, branches5, nothing,  100.0)

m=JuMP.Model()
devices_netinjection =  JumpAffineExpressionArray(length(sys5b.buses), sys5b.time_periods)

#Thermal Generator Models
p_th, inyection_array = PowerSimulations.activepowervariables(m, devices_netinjection,   sys5b.generators.thermal, sys5b.time_periods);
m = PowerSimulations.powerconstraints(m, p_th, sys5b.generators.thermal, sys5b.time_periods)
pre_set = [d for d in sys5b.generators.renewable if !isa(d, RenewableFix)]
p_re, inyection_array = PowerSimulations.activepowervariables(m, devices_netinjection,  pre_set, sys5b.time_periods)
m = PowerSimulations.powerconstraints(m, p_re, pre_set, sys5b.time_periods)
test_cl = [d for d in sys5b.loads if !isa(d, PowerSystems.StaticLoad)] # Filter StaticLoads Out
pcl, inyection_array = PowerSimulations.loadvariables(m, devices_netinjection,  test_cl, sys5b.time_periods);
m = PowerSimulations.powerconstraints(m, pcl, test_cl, sys5b.time_periods)

#Injection Array
TsNets = PowerSimulations.timeseries_netinjection(sys5b)
#CopperPlate Network test
m = PowerSimulations.copperplatebalance(m, inyection_array, TsNets, sys5b.time_periods);

#Cost Components
tl = PowerSimulations.variablecost(m, pcl, test_cl);
tre = PowerSimulations.variablecost(m, p_re, pre_set)
tth = PowerSimulations.variablecost(m, p_th, sys5b.generators.thermal);

#objective
@objective(m, Min, tl+tre+tth);


text:
game=newgame()
unit=newunit((0,0,1),(3,1,2))
placeunit!(game,unit)
@test isa(game.map[(3,1,2)],Unit)


text:
### A Pluto.jl notebook ###
# v0.14.4

using Markdown
using InteractiveUtils

# ╔═╡ 09d8b052-ac4a-11eb-09db-b3560d91e761
using Omega

# ╔═╡ 7c492773-17bf-4dd1-9b45-2ce200441bc2
using Distributions

# ╔═╡ bdd33342-30da-4f57-ba56-ec71009ee082
md"# Causal VAE"

# ╔═╡ 300a118e-cd19-4515-a103-c9623554da8d
z = 1 ~ Normal(0, 1)

# ╔═╡ 4f788865-83dc-44fc-bd60-93e129cffd28
f(θ, z) = (-θ + z, θ + z)

# ╔═╡ 2a926551-e019-4385-b60b-55682344e11f
θs = rand(5)

# ╔═╡ 9f4f2521-03ef-4baa-8588-44555c938a6e

# ╔═╡ dbcca834-f021-4c18-beaa-f6afcd172e99
function f(ω)
	z = 1 ~ Normal(0, 1)
	μ_X, σ_X = f(θs[1], z)
	x = ~ Normal(μ_X, σ_X)

	# t = 0
	μ_y_0, σ_y_0 = h(θs[2], z)
	y_0 ~ normal(mu_y_0, sigma_y_0)

	# = 1
	μ_y_1, σ_y_1 = h(θs[3], z)
	y_1 = 3 ~ Normal(μ_y_1, σ_y_1)
	return y_1 - y_0
end

# ╔═╡ 3a4faf28-565e-41b4-b278-bbb5b32c48ad
μ_X, σ_X = f(θs[1], z)

# ╔═╡ 7b60e041-3c72-44b8-ab88-4d1af04d2a8a

# ╔═╡ Cell order:
# ╟─bdd33342-30da-4f57-ba56-ec71009ee082
# ╠═09d8b052-ac4a-11eb-09db-b3560d91e761
# ╠═7c492773-17bf-4dd1-9b45-2ce200441bc2
# ╠═300a118e-cd19-4515-a103-c9623554da8d
# ╠═4f788865-83dc-44fc-bd60-93e129cffd28
# ╠═2a926551-e019-4385-b60b-55682344e11f
# ╠═3a4faf28-565e-41b4-b278-bbb5b32c48ad
# ╠═9f4f2521-03ef-4baa-8588-44555c938a6e
# ╠═dbcca834-f021-4c18-beaa-f6afcd172e99
# ╠═7b60e041-3c72-44b8-ab88-4d1af04d2a8a

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using GenericInstruments
using Test

# hardware testing done in lab. check examples folder
@testset "GenericInstruments.jl" begin
    # Write your own tests here.
end


text:
using IterableTables
using NamedTuples
using TimeSeries
using DataFrames
using NullableArrays
import TypedTables
using Base.Test

@testset "TimeSeries" begin

dates  = collect(Date(1999,1,1):Date(1999,1,3))

source_ta1 = TimeArray(dates, collect(1:length(dates)))

@test isiterable(source_ta1) == true

source_it = getiterator(source_ta1)

@test length(source_it) == 3

df1 = DataFrame(source_ta1)
@test size(df1) == (3,2)
@test df1[:timestamp] == dates
@test df1[:value] == [1,2,3]

source_ta2 = TimeArray(dates, collect(1:length(dates)), ["a"])
df2 = DataFrame(source_ta2)
@test size(df2) == (3,2)
@test df2[:timestamp] == dates
@test df2[:a] == [1,2,3]

source_ta3 = TimeArray(dates, hcat(collect(1:length(dates)), collect(length(dates):-1:1)), ["a", "b"])
df3 = DataFrame(source_ta3)
@test size(df3) == (3,3)
@test df3[:timestamp] == dates
@test df3[:a] == [1,2,3]
@test df3[:b] == [3,2,1]

source_ta4 = TimeArray(dates, hcat(collect(1:length(dates)), collect(length(dates):-1:1)))
df4 = DataFrame(source_ta4)
@test size(df4) == (3,3)
@test df4[:timestamp] == dates
@test df4[:_1] == [1,2,3]
@test df4[:_2] == [3,2,1]

source_tt = TypedTables.@Table(a=[4.,5.], time=[Date(1999,1,1),Date(1999,1,2)], b=[6.,8.], c=[12.,24.])
ta1 = TimeArray(source_tt, timestamp_column=:time)
@test length(ta1) == 2
@test TimeSeries.colnames(ta1) == ["a", "b", "c"]
@test ta1.values == [4. 6. 12.;5. 8. 24]
@test ta1.timestamp == [Date(1999,1,1),Date(1999,1,2)]

source_df = DataFrame(a=[4.,5.], timestamp=[Date(1999,1,1),Date(1999,1,2)], b=[6.,8.])
ta2 = TimeArray(source_df)
@test length(ta2) == 2
@test TimeSeries.colnames(ta2) == ["a", "b"]
@test ta2.values == [4. 6.;5. 8.]
@test ta2.timestamp == [Date(1999,1,1),Date(1999,1,2)]

end


text:
"""
    is_invariant(class::Set{String}, group::Bijection{String, T}, invariant) -> (is_invariant::Bool, values::Vector)

Verifies that the given function `invariant` is an invariant for a given (conjugacy) class.
The `class` is given by labels that match the corresponding elements in the Bijection `group`.
"""
function is_invariant(class::Set{String}, group::Bijection{String, T}, invariant; prnt=true)::Tuple{Bool, Vector} where T
    if prnt println("The calculated values of the function '$invariant' are: ") end
    values = apply_invariant_to_all_in_class(class, group, invariant; prnt)
    cond = false
    if length(values) == 1
        cond = true
    end
    return (cond, values)
end

"""
    apply_invariant_to_all_in_class(class::Set{String}, group::Bijection{String, T}, invariant; prnt=true) -> values::Vector

Calculates the values of `invariant` for every element of a given `class`.
The `class` is given by labels that match the corresponding elements in the Bijection `group`.
"""
function apply_invariant_to_all_in_class(class::Set{String}, group::Bijection{String, T}, invariant; prnt=true)::Vector where T
    values = []
    for k in class
        v = invariant(group[k])
        if !(v in values)
            if prnt print(v, ", ") end
            push!(values, v)
        end
    end
    return values
end

"""
    apply_invariant_to_first_in_class(class::Set{String}, group::Bijection{String, T}, invariant) -> value

Calculates the value of `invariant` for the first element of a given `class`.
The `class` is given by labels that match the corresponding elements in the Bijection `group`.
"""
function apply_invariant_to_first_in_class(class::Set{String}, group::Bijection{String, T}, invariant) where T
    k = first(class)
    value = invariant(group[k])
    return value
end

"""
    apply_invariant_to_first_in_all_classes(conjugacy_classes::Vector{Set{String}}, group::Bijection{String, T}, invariant; prnt=true) -> values::Vector

Calculates the values of `invariant` for the first element of all given *classes*.
The `conjugacy_classes` are given by Sets of labels that match the corresponding elements in the Bijection `group`.
"""
function apply_invariant_to_first_in_all_classes(conjugacy_classes::Vector{Set{String}}, group::Bijection{String, T}, invariant; prnt=true)::Vector where T
    values = []
    for i in 1:length(conjugacy_classes)
        class = conjugacy_classes[i]
        v = apply_invariant_to_first_in_class(class, group, invariant)
        if prnt println("Class #$i: \t$invariant(g) = $v") end
        push!(values, v)
    end
    return values
end


text:
module Caffe

using ProtoBuf
import ..Merlin
import ..Merlin.Var

include("proto/caffe_pb.jl")

function gksp(p)
  ksize = p.kernel_h > 0 ? (p.kernel_h, p.kernel_w) : (p.kernel_size, p.kernel_size)
  stride = p.stride_h > 0 ? (p.stride_h, p.stride_w) : (p.stride, p.stride)
  pad = p.pad_h > 0 ? (p.pad_h, p.pad_w) : (p.pad, p.pad)
  map(Int,ksize), map(Int,stride), map(Int,pad)
end

function conv(layer)
  blobs = layer.blobs
  w = blobs[1].data
  b = length(blobs) > 1 ? blobs[2].data : nothing
  num = blobs[1].num > 0 ? blobs[1].num : blobs[1].shape.dim[1]
  channels = blobs[1].channels > 0 ? blobs[1].channels : blobs[1].shape.dim[2]
  param = layer.convolution_param
  ksize, stride, pad = gksp(param)

  if param.group == 1
    w = reshape(w, (ksize[1],ksize[2],Int64(channels),Int64(num)))
  end
  Merlin.Conv(Merlin.Var(w), Merlin.Var(b), stride, pad)
end

function pooling(layer)
  param = layer.pooling_param
  ksize,stride,pad = gksp(param)

  if param.pool == __enum_PoolingParameter_PoolMethod().MAX
    mode = "max"
  elseif param.pool == __enum_PoolingParameter_PoolMethod().AVE
    mode = "average"
  else
    mode = nothing
  end
  Merlin.Pooling{length(pad)}(mode,ksize,stride,pad)
end

type _data
  backend
  batch_size
  crop_size
  force_encoded_color
  mean_file
  mirror
  scale
  source
end

function data(layer)
  param = layer.data_param
  _data(
    param.backend,
    param.batch_size,
    param.crop_size,
    param.force_encoded_color,
    param.mean_file,
    param.mirror,
    param.scale,
    param.source
  )
end

type _softmax_loss
  axis
end

function softmax_loss(layer)
  if isdefined(layer,:softmax_param)
    axis = layer.softmax_param.axis
  else
    axis= __val_SoftmaxParameter[:axis]
  end
  _softmax_loss(axis)
end

type _relu end

"""
Load Caffe model.
"""
function load(path)
  np = open(path) do io
    readproto(io, NetParameter())
  end
  ltype = __enum_V1LayerParameter_LayerType()
  x = Merlin.Var(:x)
  d=Any[]
  for l in (length(np.layer) > 0 ? np.layer : np.layers)
    f = nothing
    if l._type == ltype.CONVOLUTION || l._type == "Convolution"
      f = conv(l)
      #x = f(x)
    elseif l._type == ltype.POOLING || l._type == "Pooling"
      f = pooling(l)
      x = f(x)
    elseif l._type == ltype.RELU || l._type == "ReLU"
      f = _relu()
      x = Merlin.relu(x)
    elseif l._type == ltype.DATA || l._type == "Data"
      f = data(l)
      #x = f(x)
    elseif l._type == ltype.SOFTMAX_LOSS || l._type == "SoftmaxWithLoss"
      f = softmax_loss(l)
      #x = f(x)
    end
    push!(d,(l.name,typeof(f),l._type))
  end
  x,d
end

end


text:
type DiscretizedSwaption{E <: Exercise} <: DiscretizedOption
  underlying::DiscretizedSwap
  exercise::E
  exerciseTimes::Vector{Float64}
  fixedPayDates::Vector{Date}
  fixedResetDates::Vector{Date}
  floatingPayDates::Vector{Date}
  floatingResetDates::Vector{Date}
  lastPayment::Float64
  common::DiscretizedAssetCommon
end

function DiscretizedSwaption{DC <: DayCount}(swaption::Swaption, referenceDate::Date, dc::DC)
  dates = copy(swaption.exercise.dates)
  fixed_coups = swaption.swap.legs[1].coupons
  floating_coups = swaption.swap.legs[2].coupons
  nominal = swaption.swap.nominal
  swapT = swaption.swap.swapT
  n = length(dates)

  exerciseTimes = zeros(n)
  # fixedPayDates = get_pay_dates(fixed_coups)
  # fixedResetDates = get_reset_dates(fixed_coups)
  # floatingPayDates = get_pay_dates(floating_coups)
  # floatingResetDates = get_reset_dates(floating_coups)
  fixedPayDates = copy(swaption.swap.args.fixedPayDates)
  fixedResetDates = copy(swaption.swap.args.fixedResetDates)
  floatingPayDates = copy(swaption.swap.args.floatingPayDates)
  floatingResetDates = copy(swaption.swap.args.floatingResetDates)

  for i = 1:n
    exerciseTimes[i] = year_fraction(dc, referenceDate, dates[i])
  end

  # Date adjustments can get time vectors out of sync
  # Here we try and collapse similar dates which could cause a mispricing
  @simd for i = 1:n
    @inbounds exerciseDate = dates[i]

    for j = 1:length(fixed_coups)
      @inbounds if within_next_week(exerciseDate, fixedPayDates[j]) && fixedResetDates[j] < referenceDate
        @inbounds fixedPayDates[j] = exerciseDate
      end

      @inbounds if within_previous_week(exerciseDate, fixedResetDates[j])
        @inbounds fixedResetDates[j] = exerciseDate
      end
    end

    for j = 1:length(floating_coups)
      @inbounds if within_previous_week(exerciseDate, floatingResetDates[j])
        @inbounds floatingResetDates[j] = exerciseDate
      end
    end
  end

  lastFixedPayment = year_fraction(dc, referenceDate, fixedPayDates[end])
  lastFloatingPayment = year_fraction(dc, referenceDate, floatingPayDates[end])

  lastPayment = max(lastFixedPayment, lastFloatingPayment)
  underlying = DiscretizedSwap(nominal, swapT, referenceDate, dc, fixedPayDates, fixedResetDates, floatingPayDates, floatingResetDates, swaption.swap.args)
  exercise = swaption.exercise

  DiscretizedSwaption(underlying, exercise, exerciseTimes, fixedPayDates, fixedResetDates, floatingPayDates, floatingResetDates, lastPayment, DiscretizedAssetCommon())
end

function mandatory_times(discretizedSwaption::DiscretizedSwaption)
  times = mandatory_times(discretizedSwaption.underlying)
  times =  times[times .>= 0.0]
  times = vcat(times, discretizedSwaption.exerciseTimes)
  return times
end

function reset!(dSwaption::DiscretizedSwaption, sz::Int)
  initialize!(dSwaption.underlying, dSwaption.common.method, dSwaption.lastPayment)
  general_reset!(dSwaption, sz)
  return dSwaption
end


text:
module PlutoStaticHTML

using Base64: base64encode
using Pluto:
    Cell,
    CellOutput,
    Notebook,
    PlutoRunner,
    ServerSession,
    SessionActions,
    WorkspaceManager,
    generate_html,
    load_notebook_nobackup,
    update_dependency_cache!,
    update_run!,
    update_save_run!

include("html.jl")
include("build.jl")

export notebook2html, run_notebook!
export parallel_build!

end # module


text:
struct CachedEoS{EoS<:EoSModel,S} <: EoSModel
    model::EoS
    cache::Dict{Symbol,Any}
    saturation_pressure_aprox::S
end

const CACHED_SAT_PURE_APROX{S} =  CachedEoS{<:Any,S}

#is 256 too low? to high?
sizehint(::CACHED_SAT_PURE_APROX{Nothing}) = 256

function saturation_pressure_aprox_strategy(model::CachedEoS)
    return model.saturation_pressure_aprox
end
Base.length(model::CachedEoS) = Base.length(model.model)

function CachedEoS(model::EoSModel;saturation_pressure_aprox = nothing)
    cache=Dict{Symbol,Any}()
    return CachedEoS(model,cache,saturation_pressure_aprox)
end

function eos(model::CachedEoS,V,T,z=SA[1.0],phase=:unknown)
    return eos(model.model,V0,T0)
end

function eos_res(model::CachedEoS,V,T,z=SA[1.0],phase=:unknown)
    return eos_res(model.model,V,T)
end
mw(model::CachedEoS) = mw(model.model)

function Base.show(io::IO,mime::MIME"text/plain",model::CachedEoS)
    println(io,"EoS Model with cached results")
    show(io,mime,model.model)
end

function Base.show(io::IO,model::CachedEoS)
    print(io,"CachedEoS(",model.model,")")
end

function lb_volume(model::CachedEoS,z=SA[1.0])
    return lb_volume(model.model)
end

function T_scale(model::CachedEoS,z=SA[1.0])
    return T_scale(model.model)
end

function p_scale(model::CachedEoS,z=SA[1.0])
    return p_scale(model.mode)
end

function crit_pure(model::CachedEoS)
    if haskey(model.cache,:crit_pure)
        return model.cache[:crit_pure]
    else
        local res
        try
            res = crit_pure(model.model)
        catch
            res = (NaN,NaN,NaN)
        end
        model.cache[:crit_pure] = res
        return res
    end
end

function split_model(model::CachedEoS)
    if haskey(model.cache,:split_model)
        return model.cache[:split_model]
    else
        primal_splitted_model = split_model(model.model)
        res = [CachedEoS(modeli,saturation_pressure_aprox = model.saturation_pressure_aprox) for modeli in primal_splitted_model]
        model.cache[:split_model] = res
        return res
    end
end

function acentric_factor(model::CachedEoS)
    if haskey(model.cache,:acentric_factor)
        return model.cache[:acentric_factor]
    else
        local res
        try
            res = acentric_factor(model.model)
        catch
            res = NaN
        end
        model.cache[:acentric_factor] = res
        return res
    end
end

#for some reason, it does not work without an overloading here
function saturation_pressure(model::CachedEoS{<:Any},T::Real)
    aprox_strategy = saturation_pressure_aprox_strategy(model)
    return saturation_pressure_cached(aprox_strategy,model,T)
end

#no caching
function saturation_pressure_cached(::Nothing,model::CachedEoS,T)
    return saturation_pressure(model.model,T)
end

"""
    saturation_pressure!(model::CachedEoS,T,x0=x0_sat_pure(model,T))
    performs a calculation of [`saturation_pressure`](@ref) and stores the results in a vector, sorted by temperature.

"""
function saturation_pressure!(model::CachedEoS,T::S,x0 = x0_sat_pure(model,T)) where S<:Real
    #empty cache, create one
    if !haskey(model.cache,:saturation_pressure)
        (T_c, p_c, V_c) = crit_pure(model)
        T_TYPE = promote_type(S,typeof(T_c))
        T_vec,P_vec,Vl_vec,Vv_vec = T_TYPE[],T_TYPE[],T_TYPE[],T_TYPE[]
        hint = sizehint(model) #gives an expected length
        sizehint!(T_vec,hint)
        sizehint!(P_vec,hint)
        sizehint!(Vl_vec,hint)
        sizehint!(Vv_vec,hint)
        push!(T_vec,T_c)
        push!(P_vec,p_c)
        push!(Vl_vec,V_c)
        push!(Vv_vec,V_c)
        model.cache[:saturation_pressure] = (T_vec,P_vec,Vl_vec,Vv_vec)
    end
    T_vec,P_vec,Vl_vec,Vv_vec = model.cache[:saturation_pressure]
    position = searchsorted(T_vec,T)
    #one and only one match, that means that T was already evaluated
    if first(position) == last(position)
        idx = first(position)
        _P0 = convert(typeof(T),P_vec[idx])
        _Vl = convert(typeof(T),Vl_vec[idx])
        _Vv = convert(typeof(T),Vv_vec[idx])
        return _P0,_Vl,_Vv
    end
    #calculate sat pure
    T_c = last(T_vec)
    if T_c < T
        insert_first = true
    end
    P0,Vl,Vv = saturation_pressure(model.model,T,x0_sat_pure(model,T))
    #we dont want to store points too close
    insert_values = true
    if length(T_vec) >= sizehint(model)
        insert_values = false
    end

    #avoid splitting if possible
    insert_first = false
    Tmin = first(T_vec)
    if Tmin > T
        insert_first = true
    end
    if insert_values
        if insert_first
        pushfirst!(T_vec,T)
        pushfirst!(P_vec,P0)
        pushfirst!(Vl_vec,Vl)
        pushfirst!(Vv_vec,Vv)
        else
        splice!(T_vec, position, [T])
        splice!(P_vec, position, [P0])
        splice!(Vl_vec, position, [Vl])
        splice!(Vv_vec, position, [Vv])
        end
    end
    return P0,Vl,Vv
end
include("sat_pure_aprox/TPInterpolation.jl")

export CachedEoS


text:
power = [15 20 37 45 50]
probab = [0.25 0.20 0.20 0.25 0.10]

# Step 0
# Compute cost function v(w, w')

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
@with_kw mutable struct LightDark2DTarget <: AbstractLD2
    min_noise::Float64      = 0.0
    min_noise_loc::Float64  = 5.0
    term_radius::Float64    = 0.05
    init_dist::Any          = SymmetricNormal2([2.0, 2.0], 0.5)
    discount::Float64       = 1.0
    count::Int              = 0
    n_rand::Int = 0
end

reward(p::LightDark2DTarget, s::Vec2) = -1.0
reward(p::LightDark2DTarget, s::Vec2, a::Vec2) = reward(p, s)
reward(p::LightDark2DTarget, s::Vec2, a::Vec2, sp::Vec2) = reward(p, s)

# actions() with an additional boolean argument returns the available steps
# actions() without a boolean returns all combinations of steps in x and y
POMDPs.actions(p::AbstractLD2, ::Bool) = [1.0, 0.5, 0.1, 0.01, 0.0];
POMDPs.actions(p::AbstractLD2) = Vec2Iter(collect(permutations(vcat(POMDPs.actions(p, true), -POMDPs.actions(p,true)), 2)))

function isterminal(p::AbstractLD2, particles::Vector{LPDMParticle{Vec2}})
    mode_particle = LPDM.mode(particles)
    return isterminal(p, mode_particle.state)
end

# NOTE: already exists
# isterminal(p::AbstractLD2, state::Vec2) = (abs(s[1]) < p.term_radius) && (abs(s[2]) < p.term_radius)

max_stars_repo_name: ebalaban/LightDarkLPDM

max_issues_repo_name: tomerarnon/LightDarkLPDM

max_forks_repo_name: tomerarnon/LightDarkLPDM


text:
export freeze!, append!, insert_rows!, insert_cols!, delete_rows!, delete_cols!

"""
    freeze!(client::GoogleSheetsClient, sheet::Sheet, rows::Int64=0, cols::Int64=0)::Dict{Any,Any}

Freeze rows and columns in a sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `sheet::Sheet`: sheet
- `rows::Int64=0`: number of rows
- `cols::Int64=0`: number of columns
"""
function freeze!(client::GoogleSheetsClient, sheet::Sheet, rows::Int64=0, cols::Int64=0)::Dict{Any,Any}
    body = Dict(
        "requests" => [
            Dict(
                "updateSheetProperties" => Dict(
                    "properties" => Dict(
                        "sheetId" => sheet.id,
                        "gridProperties" => Dict(
                            "frozenRowCount" => rows,
                        ),
                    ),
                    "fields" => "gridProperties.frozenRowCount",
                ),
            ),

            Dict(
                "updateSheetProperties" => Dict(
                    "properties" => Dict(
                        "sheetId" => sheet.id,
                        "gridProperties" => Dict(
                            "frozenColumnCount" => cols,
                        ),
                    ),
                    "fields" => "gridProperties.frozenColumnCount",
                ),
            ),
        ],
    )

    return batch_update!(client, sheet.spreadsheet, body)
end

"""
    append!(client::GoogleSheetsClient, sheet::Sheet, rows::Int64=0, cols::Int64=0)::Dict{Any,Any}

Append rows and columns to a sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `sheet::Sheet`: sheet
- `rows::Int64=0`: number of rows
- `cols::Int64=0`: number of columns
"""
function Base.append!(client::GoogleSheetsClient, sheet::Sheet, rows::Int64=0, cols::Int64=0)::Dict{Any,Any}
    body = Dict(
        "requests" => [
            Dict(
                "appendDimension" => Dict(
                    "sheetId" => sheet.id,
                    "dimension" => "ROWS",
                    "length" => rows,
                ),
            ),

            Dict(
                "appendDimension" => Dict(
                    "sheetId" => sheet.id,
                    "dimension" => "COLUMNS",
                    "length" => cols,
                ),
            ),
        ],
    )

    return batch_update!(client, sheet.spreadsheet, body)
end

"""
    insert_rows!(client::GoogleSheetsClient, range::CellIndexRange1D)::Dict{Any,Any}

Insert rows into to a sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `range::CellIndexRange1D`: cell index range
"""
function insert_rows!(client::GoogleSheetsClient, range::CellIndexRange1D)::Dict{Any,Any}
    return _insert!(client, range, "ROWS")
end

"""
    insert_cols!(client::GoogleSheetsClient, range::CellIndexRange1D)::Dict{Any,Any}

Insert columns into to a sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `range::CellIndexRange1D`: cell index range
"""
function insert_cols!(client::GoogleSheetsClient, range::CellIndexRange1D)::Dict{Any,Any}
    return _insert!(client, range, "COLUMNS")
end

"""
Insert rows or columns into to a sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `range::CellIndexRange1D`: cell index range
- `dim::AbstractString`: dimension string
"""
function _insert!(client::GoogleSheetsClient, range::CellIndexRange1D, dim::AbstractString)::Dict{Any,Any}
    body = Dict(
        "requests" => [
            Dict(
                "insertDimension" => Dict(
                    "range" => gsheet_json(range, dim)
                ),
            ),
        ],
    )

    return batch_update!(client, range.sheet.spreadsheet, body)
end

"""
Delete rows from a sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `range::CellIndexRange1D`: cell index range
"""
function delete_rows!(client::GoogleSheetsClient, range::CellIndexRange1D)::Dict{Any,Any}
    return _delete!(client, range, "ROWS")
end

"""
    delete_cols!(client::GoogleSheetsClient, range::CellIndexRange1D)::Dict{Any,Any}

Delete columns from a sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `range::CellIndexRange1D`: cell index range
"""
function delete_cols!(client::GoogleSheetsClient, range::CellIndexRange1D)::Dict{Any,Any}
    return _delete!(client, range, "COLUMNS")
end

"""
Delete rows or columns from a sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `range::CellIndexRange1D`: cell index range
- `dim::AbstractString`: dimension string
"""
function _delete!(client::GoogleSheetsClient, range::CellIndexRange1D, dim::AbstractString)::Dict{Any,Any}
    body = Dict(
        "requests" => [
            Dict(
                "deleteDimension" => Dict(
                    "range" => gsheet_json(range, dim)
                ),
            ),
        ],
    )

    return batch_update!(client, range.sheet.spreadsheet, body)
end


text:
using JSON
using DataStructures: Deque
#Taps = include("/home/schinavro/libCalc/taps/taps_parallel/taps.jl")
#using .Taps: Coords

int, float, real, number, bool = 1, 2, 3, 4, 5
int8, int16, int32, int64 = 6, 7, 8, 9
uint8, uint16, uint32, uint64 = 10, 11, 12, 13
float16, float32, float64 = 14, 15, 16
complex64, complex128 = 17, 18

const Float = AbstractFloat

typemaps = Dict(
    int=>Int, Int=>int,
    float=>Float, Float=>float, Real=>float, Number=>float,
    bool=>Bool, Bool=>bool,
    int8=>Int8, Int8=>int8,
    int16=>Int16, Int16=>int16,
    int32=>Int32, Int32=>int32,
    int64=>Int64, Int64=>int64,
    uint8=>UInt8, UInt8=>uint8,
    uint16=>UInt16, UInt16=>uint16,
    uint32=>UInt32, UInt32=>uint32,
    uint64=>UInt64, UInt64=>uint64,
    float16=>Float16, Float16=>float16,
    float32=>Float32, Float32=>float32,
    float64=>Float64, Float64=>float64,
    complex64=>Complex{Float32}, Complex{Float32}=>complex64,
    complex128=>Complex{Float64}, Complex{Float64}=>complex128
)
ordermaps = Dict("C"=>0, 0=>"C", 1=>"F", "F"=>1)

function write_header(pointer, arr)
    header = []
    dtype = eltype(arr)
    ndim = ndims(arr)
    shape = size(arr)
    order = "F"
    header = Array{Int64, 1}([pointer, typemaps[dtype], ndim, shape..., ordermaps[order]])
    headerbytes = reinterpret(UInt8, header)
    return headerbytes
end

function read_header(arrbytes)
    header_size = reinterpret(Int64, arrbytes[1:8])[1]
    header = reinterpret(Int64, arrbytes[9:8+header_size])
    pointer, dtype, ndim = header[1], typemaps[header[2]], header[3]
    shape = header[4:3+ndim]
    order = ordermaps[header[4+ndim]]
    header = header_size, pointer, dtype, ndim, shape, order
    return header
end

function statify(arrlist, d::Union{Dict, Array{Any, 1}})
    pointerlist = ["__$(arr[1])__"  for arr in arrlist]
    global queue = Deque{Tuple}()
    push!(queue, (hash(d), d))
    memo = Set()
    while true
        global queue
        isempty(queue) ? break : nothing
        id_, o = popfirst!(queue)
        id_ in memo ? continue : nothing
        union!(memo, id_)

        if typeof(o)<:Dict for (k, v) in pairs(o)
            push!(queue, (hash(v), v))
            v in pointerlist ? o[k] = arrlist[findfirst(x->x==v, pointerlist)][2] : nothing
            end
        elseif typeof(o)<:Union{Array{Any, 1}, Tuple} for i=1:length(o)
            v = o[i]
            push!(queue, (hash(v), v))
            v in pointerlist ? o[i] = arrlist[findfirst(x->x==v, pointerlist)][2] : nothing
            end
        end
    end
    return d
end

function pointify(d, pointer, binarylist::Array)
    global queue = Deque{Tuple}()
    push!(queue, (hash(d), d))
    memo = Set()
    while true
        global queue
        isempty(queue) ? break : nothing
        id_, o = popfirst!(queue)
        id_ in memo ? continue : nothing
        union!(memo, id_)
        if typeof(o)<:Dict for (k, v) in pairs(o)
            push!(queue, ((hash(v), v)))
            if typeof(v)<:Array{T, N} where {T<:Number, N}
                pstr = "__$(pointer)__"

                headerbytes = write_header(pointer, o[k])
                binar = cat(reinterpret(UInt8, [length(headerbytes)]), headerbytes, reinterpret(UInt8, vec(v)), dims=1)
                push!(binarylist, binar)
                o[k] = pstr
                pointer += 1
            end end
        elseif typeof(o)<:Union{Array{Any, 1}, Tuple} for i=1:length(o)
            v = o[i]
            push!(queue, (hash(v), v))
            if typeof(v)<:Array{T, N} where {T<:Number, N}

                pstr = "__$(pointer)__"
                headerbytes = write_header(pointer, v)
                binar = cat(reinterpret(UInt8, [length(headerbytes)]), headerbytes, reinterpret(UInt8, vec(v)), dims=1)
                push!(binarylist, binar)
                o[i] = pstr
                pointer += 1
            end end
        end
    end
    return d, pointer, binarylist
end

function packing(args...; kwargs...)::Array{UInt8, 1}
    args, kwargs = Array{Any, 1}([args...]), Dict{Symbol, Any}(kwargs)
    binarylist = []
    pointer = 0

    args, pointer, binarylist = pointify(args, pointer, binarylist)
    kwargs, pointer, binarylist = pointify(kwargs, pointer, binarylist)
    kwargs = Dict{Any, Any}(kwargs)

    kwargs[:args] = args

    howmanybinary = binarylist == [] ? 0 : reinterpret(UInt8, [length(binarylist)])

    eachsize = binarylist == [] ? [] : reinterpret(UInt8, [length(b) for b in binarylist])

    binarybytes = cat(howmanybinary, eachsize, cat(binarylist..., dims=1), dims=1)
    kwargsbytes = Array{UInt8}(JSON.json(kwargs))

    data = cat(binarybytes, kwargsbytes, dims=1)
    header = reinterpret(UInt8, [length(data)+16, length(binarybytes), length(kwargsbytes)])
    return cat(header, data, dims=1)
end

function unpacking(bytesarr::Array{UInt8, 1})
    binarysize, kwargssize = reinterpret(Int64, bytesarr[1:16])
    partition = 16 + binarysize

    binarybytes = bytesarr[17:partition]
    kwargsbytes = bytesarr[partition+1:end]

    #binarybytes to binarylist
    howmanybinary = reinterpret(Int64, binarybytes[1:8])[1]
    eachsize = reinterpret(Int64, binarybytes[9:8 + 8howmanybinary])
    partition = 9 + 8howmanybinary

    arrlist = []
    for size in eachsize
        arrbytes = binarybytes[partition:partition+size-1]
        header_size, pointer, dtype, ndim, shape, order = read_header(arrbytes)

        shape = order == "C" ? reverse(shape) : shape

        arr = reinterpret(dtype, arrbytes[header_size+9:end])

        array = Array{dtype, length(shape)}(reshape(arr, shape...))
        push!(arrlist, (pointer, array))
        partition += size
    end

    pkwargs = JSON.parse(String(kwargsbytes); dicttype=Dict{Symbol,Any})
    pargs = pkwargs[:args]
    delete!(pkwargs, :args)

    # link binarybytes and everything
    args = statify(arrlist, pargs)
    kwargs = statify(arrlist, pkwargs)

    return args, kwargs
end

max_stars_repo_name: schinavro/taps

max_issues_repo_name: schinavro/taps

max_forks_repo_name: schinavro/taps


text:
facts("Bounding Box") do
  context("box contains point") do
    @fact contains(box1, Point(5,5)) => true
    @fact contains(box1, Point(-5,5)) => false
    @fact contains(box1, Point(5,15)) => false
    @fact contains(box1, Point(10,9)) => false
    @fact contains(box1, Point(11,11)) => false
  end

  context("box intersects box") do
    @fact intersects(box1, box1) => true
    @fact intersects(box1, box2) => true
    @fact intersects(box1, box3) => true
    @fact intersects(box1, box4) => false
    @fact intersects(box1, box5) => false
    @fact intersects(box1, box6) => true
  end

  context("box contains box") do
    @fact contains(box1, box1) => false
    @fact contains(box1, box2) => true
    @fact contains(box1, box3) => false
    @fact contains(box1, box4) => false
    @fact contains(box1, box5) => false
    @fact contains(box1, box6) => false
  end

  context("box intersect box") do
    @fact isequal(intersect(box1, box1), Nullable(box1)) => true
    @fact isequal(intersect(box1, box2), Nullable(box2)) => true
    @fact isequal(intersect(box1, box3), Nullable(Box{Int}(Point(1,1), Point(10,10)))) => true
    @fact isequal(intersect(box1, box4), Nullable{Box{Int}}()) => true
    @fact isequal(intersect(box1, box5), Nullable{Box{Int}}()) => true
    @fact isequal(intersect(box1, box6), Nullable(box1)) => true
  end
end


text:
@testset "Simplify" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  K, a = NumberField(x^2 - 10, cached = false)
  L, mL = simplify(K)
  @test is_isomorphic(K, L)

  K, a = NumberField(x^2 - 100000000000000000000000, cached = false)
  L, mL = simplify(K)
  @test is_isomorphic(K, L)

  K, a = NumberField(x - 10, cached = false)
  L, mL = simplify(K)
  @test is_isomorphic(K, L)

  K,a = NumberField(x^4-100020001*x^2+100040006000400010, cached = false)
  L, mL = simplify(K)
  @test is_isomorphic(K, L)

  K, a = number_field(8*x^3 + 4*x^2 - 4*x - 1, cached = false)
  L, mL = simplify(K)
  @test Hecke.is_defining_polynomial_nice(L)
  L1, mL1 = simplify(K, canonical = true)
  @test Hecke.is_defining_polynomial_nice(L1)
  @test is_isomorphic(L1, L)

  f = -1//3*x^2 - x + 1
  K, a = number_field(f, cached = false)
  L, mL = simplify(K)
  @test is_isomorphic(K, L)

  f = x - 1
  K, a = number_field(f, cached = false)
  _, g = Hecke.polredabs(K)
  @test g == f
end


text:
export IMTFile, IMTFileCalib, IMTFileMeas, addTrailingSingleton

abstract type IMTFile <: MPIFile end

# We use a dedicated type for calib and meas. If both versions
# are the same we use the abstract type IMTFile
mutable struct IMTFileCalib <: IMTFile
  filename::String
  file::HDF5File
  mmap_measData
end

#IMTFileCalib(filename::String, file=h5open(filename,"r")) =
#   IMTFileCalib(filename, file, nothing)

function IMTFileCalib(filename::String, file=h5open(filename,"r"))
  f = IMTFileCalib(filename, file, nothing)

 return f
end

mutable struct IMTFileMeas <: IMTFile
  filename::String
  file::HDF5File
  mmap_measData
end

function IMTFileMeas(filename::String, file=h5open(filename,"r"))
  f = IMTFileMeas(filename, file, nothing)

  return f
end

# Dispatch on file extension
function (::Type{IMTFile})(filename::String, file = h5open(filename,"r"))
 if !exists(file, "/measurements")
   return IMTFileCalib(filename, file)
 else
   return IMTFileMeas(filename, file)
 end
end

function getindex(f::IMTFile, parameter)
  if exists(f.file, parameter)
    return read(f.file, parameter)
  else
    return nothing
  end
end

function getindex(f::IMTFile, parameter, default)
  if exists(f.file, parameter)
    return read(f.file, parameter)
  else
    return default
  end
end

# general parameters
version(f::IMTFile) = v"0.0.0"
uuid(f::IMTFile) = uuid4()
time(f::IMTFile) = Dates.unix2datetime(0)

# study parameters
studyName(f::IMTFile) = "n.a."
studyNumber(f::IMTFile) = 0
studyUuid(f::IMTFile) = uuid4()
studyDescription(f::IMTFile) = "n.a."
studyTime(f::IMTFile) = nothing

# experiment parameters
experimentName(f::IMTFile) = "n.a."
experimentNumber(f::IMTFile) = 0
experimentUuid(f::IMTFile) = uuid4()
experimentDescription(f::IMTFile) = "n.a."
experimentSubject(f::IMTFile) = "n.a."
experimentIsSimulation(f::IMTFile) = true
experimentIsCalibration(f::IMTFileMeas) = false
experimentIsCalibration(f::IMTFileCalib) = true
experimentHasReconstruction(f::IMTFile) = false
experimentHasMeasurement(f::IMTFile) = true

# tracer parameters
tracerName(f::IMTFile)::Vector{String} = _makeStringArray(["n.a."])
tracerBatch(f::IMTFile)::Vector{String} = ["n.a."]
tracerVolume(f::IMTFile)::Vector{Float64} = [0.0]
tracerConcentration(f::IMTFile)::Vector{Float64} = [0.0]
tracerSolute(f::IMTFile)::Vector{String} = ["Fe"]
tracerInjectionTime(f::IMTFile) = [Dates.unix2datetime(0)]
tracerVendor(f::IMTFile)::Vector{String} = ["n.a."]

# scanner parameters
scannerFacility(f::IMTFile)::String = "n.a."
scannerOperator(f::IMTFile)::String = "n.a."
scannerManufacturer(f::IMTFile)::String = "n.a."
scannerName(f::IMTFile)::String = "n.a."
scannerTopology(f::IMTFile)::String = "n.a."

# acquisition parameters
acqStartTime(f::IMTFile)::DateTime = Dates.unix2datetime(0)
acqNumAverages(f::IMTFileCalib)::Int = 1
acqNumAverages(f::IMTFileMeas)::Int = 1
acqNumFrames(f::IMTFileCalib)::Int = 1
acqNumFrames(f::IMTFileMeas)::Int = 1
acqNumPeriodsPerFrame(f::IMTFile)::Int = 1

acqGradient(f::IMTFile)::Array{Float64,4} = reshape(Matrix(Diagonal([0.0,0.0,0.0])), 3,3,1,1)
acqOffsetField(f::IMTFile)::Array{Float64,3} = reshape([0.0,0.0,0.0],3,1,1)

# drive-field parameters
dfNumChannels(f::IMTFileMeas) = size(f["/measurements"], 2)
dfNumChannels(f::IMTFileCalib) = size(f["/numberOfAvailableFrequencies"],1)
dfStrength(f::IMTFile) = [0.0 0.0 0.0] # addTrailingSingleton( addLeadingSingleton(f["/acquisition/drivefield/strength"], 2), 3)
dfPhase(f::IMTFile) = [0.0 0.0 0.0]
dfBaseFrequency(f::IMTFile) = 2.5e6
dfCustomWaveform(f::IMTFile) = "n.a."
dfDivider(f::IMTFile) = reshape([102; 96; 99],:,1)
dfWaveform(f::IMTFile) = "sine"
dfCycle(f::IMTFile) = f["/timeLength"][1]

# receiver parameters
rxNumChannels(f::IMTFileMeas) = size(f["/measurements"],2)
rxNumChannels(f::IMTFileCalib) = length(f["/numberOfAvailableFrequencies"])
rxBandwidth(f::IMTFile)::Float64 = 1.25e6
rxNumSamplingPoints(f::IMTFile) = (f["/numberOfAvailableFrequencies"][1]-1)*2
rxTransferFunction(f::IMTFile) = nothing
rxInductionFactor(f::IMTFile) = nothing

rxUnit(f::IMTFile) = "a.u."
rxDataConversionFactor(f::IMTFile) = repeat([1.0, 0.0], outer=(1,rxNumChannels(f)))
#rxDataConversionFactor(f::IMTFileMeas) = f["/acquisition/receiver/dataConversionFactor"]

# measurements
function measData(f::IMTFile, frames=1:acqNumFrames(f), periods=1:acqNumPeriodsPerFrame(f),
                  receivers=1:rxNumChannels(f))

  if !exists(f.file, "/measurements")
    # file is calibration
    dataFD = f["/systemResponseFrequencies"]
    dataFD = reshape(reinterpret(Complex{eltype(dataFD)}, vec(dataFD)), (div(size(dataFD,1),2),size(dataFD,2),size(dataFD,3),size(dataFD,4)))
    return dataFD = reshape(dataFD, (size(dataFD,1)*size(dataFD,2)*size(dataFD,3), div(size(dataFD,4),2), length(receivers), length(frames)))
  end

  tdExists = exists(f.file, "/measurements")

  if tdExists
    dataTD = f["/measurements"]
    #TODO implement for frames > 1
    ##dataFD = rfft(reshape(dataTD, size(dataTD,1), size(dataTD,2), 1, length(frames)))
    dataTD = reshape(dataTD, size(dataTD,1), size(dataTD,2), 1, length(frames))
    ##dataFD = reshape(dataTD, size(dataTD,1), size(dataTD,2), length(frames))
    return dataTD
  end
end

function systemMatrix(f::IMTFileCalib, rows, bgCorrection=true)

  if !experimentIsCalibration(f)
    return nothing
  end
  if f.mmap_measData == nothing
    f.mmap_measData = readmmap(f.file["/systemResponseFrequencies"])
  end

  data = f.mmap_measData[:, :, :, rows]
  return reshape(reinterpret(Complex{eltype(data)}, vec(data)), :, size(data,4))
end

measIsFourierTransformed(f::IMTFileMeas) = false
measIsFourierTransformed(f::IMTFileCalib) = true

measIsTFCorrected(f::IMTFile) = false
measIsSpectralLeakageCorrected(f::IMTFile) = false

measIsBGCorrected(f::IMTFile) = false

measIsFrequencySelection(f::IMTFile) = false

measIsFastFrameAxis(f::IMTFileCalib) = true
measIsFastFrameAxis(f::IMTFileMeas) = false

measIsFramePermutation(f::IMTFileCalib) = false
measIsFramePermutation(f::IMTFileMeas) = false

measIsBGFrame(f::IMTFile) = zeros(Bool, acqNumFrames(f))

measFramePermutation(f::IMTFileCalib) = nothing
measFramePermutation(f::IMTFileMeas) = nothing

#fullFramePermutation(f::IMTFile) = fullFramePermutation(f, calibIsMeanderingGrid(f))

#calibrations
calibSNR(f::IMTFileCalib) = 100*ones(rxNumFrequencies(f),rxNumChannels(f),1)
calibFov(f::IMTFile) = f["/fov"]
calibFovCenter(f::IMTFile) = [0.0,0.0,0.0]
calibSize(f::IMTFile) = [div(size(f["/systemResponseFrequencies"],1),2),
                         size(f["/systemResponseFrequencies"],2),
                         size(f["/systemResponseFrequencies"],3)]
calibOrder(f::IMTFile) = "xyz"
calibOffsetField(f::IMTFileCalib) = nothing
calibDeltaSampleSize(f::IMTFile) = [0.0,0.0,0.0]
calibMethod(f::IMTFile) = "simulation"
#calibIsMeanderingGrid(f::IMTFile) = Bool(f["/calibration/isMeanderingGrid", 0])
calibPositions(f::IMTFileCalib) = f["/calibration/positions"]

# additional functions that should be implemented by an MPIFile
filepath(f::IMTFile) = f.filename


text:
export BrickletMotionDetectorIdentity
struct BrickletMotionDetectorIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletMotionDetector
"""
Passive infrared (PIR) motion sensor with 7m range
"""
mutable struct BrickletMotionDetector <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletMotionDetector(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_motion_detector")
        deviceInternal = package.BrickletMotionDetector(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_motion_detected
"""
    $(SIGNATURES)

Returns 1 if a motion was detected. How long this returns 1 after a motion
was detected can be adjusted with one of the small potentiometers on the
Motion Detector Bricklet, see :ref:`here
<motion_detector_bricklet_sensitivity_delay_block_time>`.

There is also a blue LED on the Bricklet that is on as long as the Bricklet is
in the "motion detected" state.
"""
function get_motion_detected(device::BrickletMotionDetector)
    return device.deviceInternal.get_motion_detected()
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status led configuration.

By default the status LED turns on if a motion is detected and off is no motion
is detected.

You can also turn the LED permanently on/off.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function set_status_led_config(device::BrickletMotionDetector, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function get_status_led_config(device::BrickletMotionDetector)
    return device.deviceInternal.get_status_led_config()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletMotionDetector)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletMotionDetector, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
# test safeguarding related functions

using COSMO, SparseArrays, LinearAlgebra, Test

@testset "Acceleration safeguarding" begin
    q = [1; 1.];
    P = sparse([4. 1; 1 2]);
    A = [1. 1; 1 0; 0 1];
    l = [1.; 0; 0];
    u = [1; 0.7; 0.7];

    Aa = [-A; A]
    ba = [u; -l]
    constraint1 = COSMO.Constraint(Aa, ba, COSMO.Nonnegatives);
    m, n = size(Aa)
    # Next, we define the settings object, the model and then assemble everything:
    # model = COSMO.Model();
    # assemble!(model, P, q, constraint1, settings = COSMO.Settings(max_iter = 10, scaling = 0));
	# model.rws = COSMO.ResidualWorkspace{Float64}(m, n, model.p.C)
    # rws = model.rws
    # # make sure rws.ν works as a view
    # @. rws.sol = rand()
    # @test rws.sol[n+1:end] == rws.ν

    # # do ten iterations
    # result = COSMO.optimize!(model);
    # w0 = deepcopy(model.vars.w)
    # # do one more iteratios
    # model.settings.max_iter = 1
    # result = COSMO.optimize!(model);
    # w1 = deepcopy(model.vars.w)

    # # now put w0 through fixed-point residual function
    # fp_res = COSMO.fixed_point_residual_norm(model.rws, model, w0)
    # @test fp_res == norm(w0 - w1, 2)

    # now check algorithm behaviour with safeguarding
    model = COSMO.Model();
    accelerator = with_options(AndersonAccelerator{Float64, NoRegularizer, Type2, RollingMemory}, mem = 5, safeguarded = true, activation_reason = ImmediateActivation())
    settings = COSMO.Settings(eps_abs = 1e-6, eps_rel = 1e-6, accelerator = accelerator)
    assemble!(model, P, q, constraint1, settings = settings);
    result = COSMO.optimize!(model);

end


text:
# setup path -
_PATH_TO_PROJECT = pwd()
_PATH_TO_SRC = joinpath("execute_me_PS2.jl", "src")
_PATH_TO_CONFIG = joinpath("execute_me_PS2.jl", "config")

# import PKG, activate the project and instaintiate -
import Pkg
Pkg.activate("execute_me_PS2.jl")
Pkg.instantiate()

# load packages -
using DifferentialEquations
using Plots
using LaTeXStrings
using TOML
using Logging

# my codes -

# Problem.jl
path_to_problem = joinpath(src,"Problem.jl")
include("Problem.jl")

# Balances.jl
path_to_balances = joinpath(src,"Balances.jl")
include("Balances.jl")

# Control.jl
path_to_control = joinpath(src,"Control.jl")
include("Control.jl")

# Kinetics.jl
path_to_kinetics = joinpath(src,"Kinetics.jl")
include("Kinetics.jl")

# Solver.jl
path_to_solver = joinpath(src,"Solver.jl")
include("Solver.jl")


text:
abstract type Spline end

struct NearestNeighborSpline <: Spline
    a::Vector{Float64}
    x::Vector{Float64}
end
struct LinearSpline <: Spline
    b::Vector{Float64}
    a::Vector{Float64}
    x::Vector{Float64}
end

struct QuadraticSpline <: Spline
    c::Vector{Float64}
    b::Vector{Float64}
    a::Vector{Float64}
    x::Vector{Float64}
end

struct CubicSpline <: Spline
    d::Vector{Float64}
    c::Vector{Float64}
    b::Vector{Float64}
    a::Vector{Float64}
    x::Vector{Float64}
end

struct HermiteSpline <: Spline
    b::Vector{Float64}
    a::Vector{Float64}
    x::Vector{Float64}
end


text:
using Fretboards
using Fretboards: parse_pitch
using Test

@testset "pitch parsing and representation" begin
    maybe0 = ["", "0"]
    sharps = ["♯", "#"]
    flats = ["♭", "b"]

    "All combinations of the first letter and vectors of strings."
    function all_combinations(letter::Char, other_options...)
        str = string(letter)
        mapfoldl(x -> x isa AbstractString ? [x] : x,
                 (x, y) -> vec(x .* permutedims(y)),
                 other_options;
                 init = [lowercase(str), uppercase(str)])
    end

    @test all(parse_pitch.(all_combinations('c', maybe0)) .== Pitch(0))
    @test repr(Pitch(0)) == "pitch\"C0\""
    @test all(parse_pitch.(all_combinations('c', sharps, maybe0)) .== Pitch(1))
    @test repr(Pitch(1)) == "pitch\"C♯0\""
    @test all(parse_pitch.(all_combinations('c', flats, maybe0)) .== Pitch(-1))
    @test repr(Pitch(-1)) == "pitch\"B-1\""
    @test all(parse_pitch.(all_combinations('A', "2")) .== Pitch(24 + 9))
    @test repr(Pitch(24 + 9)) == "pitch\"A2\""
    @test all(parse_pitch.(all_combinations('A', flats, "2")) .== Pitch(24 + 8))
    @test repr(Pitch(24 + 8)) == "pitch\"G♯2\""
    @test all(parse_pitch.(all_combinations('A', sharps, "2")) .== Pitch(24 + 10))
    @test repr(Pitch(24 + 10)) == "pitch\"A♯2\""

    @test_throws ArgumentError parse_pitch("x")    # invalid pitch
    @test_throws ArgumentError parse_pitch("ax")   # trailing clutter
    @test_throws ArgumentError parse_pitch("a-1x") # trailing clutter
end

@testset "pitch calculations" begin
    @test Pitch(20) - Pitch(10) == 10
    @test Pitch(20) - Pitch(10) == 10
    @test Pitch(10) + 20 == Pitch(30)
    @test Pitch(10) - 20 == Pitch(-10)
    @test pitch"A♯" ≂ pitch"A♯-1" ≂ pitch"A♯2"
    @test !(pitch"C" ≂ pitch"D")
end


text:
abstract DummyTarget

kernel_queue = []

dummy_call(f, args) = f(args...)

function translate_kernel(target::Type{DummyTarget}, kernel_id, signature)
    info("Translating $kernel_id $signature")

    # Get the typed AST to translate:
    functype, argtypes = signature
    ast = code_typed(functype, argtypes...)

    info("Unboxing variables...")
    ast = unbox(ast)

    info("Recovering structured control flow...")
    flow = FlowGraph(ast)
    showgraph(flow)

    statements, _ = raise_flow(1, flow)
    for s in statements
        println(s)
    end

    info("Done translating.")

    dummy_call
end


text:
"""
`PA_3197`
# Frequency (MHz)
# Measured Sensitivity (mV/MPa)
# Measured Uncertainty (%)
"""

@calibration :PA_3197 [
    1 34 18
    2 51 18
    3 59 18
    4 60 18
    5 53 18
    6 52 18
    7 58 17
    8 45 18
    9 43 18
    10 43 17
    11 46 18
    12 46 18
    13 48 18
    14 49 18
    15 49 18
    16 49 18
    17 48 19
    18 48 20
    19 47 19
    20 47 20
    21 48 20
    22 48 20
    23 50 20
    24 49 20
    25 46 20
    26 47 20
    27 44 20
    28 45 20
    29 43 21
    30 43 20
]


text:
module RPrompt

using REPL
import REPL: REPL, LineEdit, REPLCompletions
import ..REvalError
import ..Const
import ..RCall:
    libR,
    rparse_p,
    reval_p,
    findNamespace,
    rcall_p,
    rprint,
    rcopy,
    render,
    protect,
    unprotect,
    prepare_inline_julia_code,
    RParseIncomplete,
    RException,
    RParseError,
    REvalError,
    RParseEOF

function simple_showerror(io::IO, er)
    Base.with_output_color(:red, io) do io
        print(io, "ERROR: ")
        showerror(io, er)
        println(io)
    end
end

function parse_status(script::String)
    status = :ok
    try
        render(script)
    catch ex
        if isa(ex, RParseIncomplete)
            status = :incomplete
        else
            status = :error
        end
    end
    status
end

function repl_eval(script::String, stdout::IO, stderr::IO)
    local nprotect = 0
    try
        script, symdict = render(script)
        if length(symdict) > 0
            Core.eval(Main, prepare_inline_julia_code(symdict))
        end
        ret = protect(reval_p(rparse_p("withVisible({$script})"), Const.GlobalEnv.p))
        nprotect += 1
        # print if the last expression is visible
        if rcopy(Bool, ret[:visible])
             rprint(stdout, ret[:value])
        end
    catch ex
        if isa(ex, REvalError)
            println(stderr, ex.msg)
        elseif isa(ex, RParseIncomplete) || isa(ex, RParseError)  || isa(ex, RParseEOF)
            println(stderr, ex.msg)
        else
            simple_showerror(stderr, ex)
        end
    finally
        unprotect(nprotect)
        return nothing
    end
end

function bracketed_paste_callback(s, o...)
    input = LineEdit.bracketed_paste(s)
    sbuffer = LineEdit.buffer(s)
    curspos = position(sbuffer)
    seek(sbuffer, 0)
    shouldeval = (bytesavailable(sbuffer) == curspos && search(sbuffer, UInt8('\n')) == 0)
    seek(sbuffer, curspos)
    if curspos == 0
        # if pasting at the beginning, strip leading whitespace
        input = lstrip(input)
    end

    if !shouldeval
        LineEdit.edit_insert(s, input)
        return
    end

    LineEdit.edit_insert(sbuffer, input)
    input = String(take!(sbuffer))

    oldpos = start(input)
    nextpos = 0
    # parse the input line by line
    while !done(input, oldpos)
        nextpos = search(input, '\n', nextpos+1)
        if nextpos == 0
            nextpos = endof(input)
        end
        block = input[oldpos:nextpos]
        status = parse_status(block)

        if status == :error  || (status == :incomplete && done(input, nextpos+1)) ||
                (done(input, nextpos+1) && !endswith(input, '\n'))
            # error / continue but not the end / at the end but no new line
            LineEdit.replace_line(s, input[oldpos:end])
            LineEdit.refresh_line(s)
            break
        elseif status == :incomplete && !done(input, nextpos+1)
            continue
        end

        if !isempty(strip(block))
            # put the line on the screen and history
            LineEdit.replace_line(s, strip(block))
            LineEdit.commit_line(s)
            # execute the statement
            terminal = LineEdit.terminal(s)
            REPL.raw!(terminal, false) && LineEdit.disable_bracketed_paste(terminal)
            LineEdit.mode(s).on_done(s, LineEdit.buffer(s), true)
            REPL.raw!(terminal, true) && LineEdit.enable_bracketed_paste(terminal)
        end
        oldpos = nextpos + 1
    end
    LineEdit.refresh_line(s)
end

mutable struct RCompletionProvider <: LineEdit.CompletionProvider
    r::REPL.LineEditREPL
end

function LineEdit.complete_line(c::RCompletionProvider, s)
    buf = s.input_buffer
    partial = String(buf.data[1:buf.ptr-1])
    # complete latex
    full = LineEdit.input_string(s)
    ret, range, should_complete = REPLCompletions.bslash_completions(full, lastindex(partial))[2]
    if length(ret) > 0 && should_complete
        return map(REPLCompletions.completion_text, ret), partial[range], should_complete
    end

    # complete r
    utils = findNamespace("utils")
    rcall_p(utils[".assignLinebuffer"], partial)
    rcall_p(utils[".assignEnd"], length(partial))
    token = rcopy(rcall_p(utils[".guessTokenFromLine"]))
    rcall_p(utils[".completeToken"])
    ret = rcopy(Array, rcall_p(utils[".retrieveCompletions"]))
    if length(ret) > 0
        return ret, token, true
    end

    return String[], 0:-1, false
end

function create_r_repl(repl, main)
    r_mode = LineEdit.Prompt("R> ";
        prompt_prefix=Base.text_colors[:blue],
        prompt_suffix=main.prompt_suffix,
        sticky=true)

    hp = main.hist
    hp.mode_mapping[:r] = r_mode
    r_mode.hist = hp
    r_mode.complete = RCompletionProvider(repl)
    r_mode.on_enter = (s) -> begin
        status = parse_status(String(take!(copy(LineEdit.buffer(s)))))
        status == :ok || status == :error
    end
    r_mode.on_done = (s, buf, ok) -> begin
        if !ok
            return REPL.transition(s, :abort)
        end
        script = String(take!(buf))
        if !isempty(strip(script))
            REPL.reset(repl)
            try
                repl_eval(script, repl.t.out_stream, repl.t.err_stream)
            catch y
                # should never reach
                simple_showerror(repl.t.err_stream, y)
            end
        end
        REPL.prepare_next(repl)
        REPL.reset_state(s)
        s.current_mode.sticky || REPL.transition(s, main)
    end

    bracketed_paste_mode_keymap = Dict{Any,Any}(
        "\e[200~" => bracketed_paste_callback
    )

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    prefix_prompt, prefix_keymap = LineEdit.setup_prefix_keymap(hp, r_mode)

    mk = REPL.mode_keymap(main)
    # ^C should not exit prompt
    delete!(mk, "^C")

    b = Dict{Any,Any}[
        bracketed_paste_mode_keymap,
        skeymap, mk, prefix_keymap, LineEdit.history_keymap,
        LineEdit.default_keymap, LineEdit.escape_defaults
    ]
    r_mode.keymap_dict = LineEdit.keymap(b)

    r_mode
end

function repl_init(repl)
    mirepl = isdefined(repl,:mi) ? repl.mi : repl
    main_mode = mirepl.interface.modes[1]
    r_mode = create_r_repl(mirepl, main_mode)
    push!(mirepl.interface.modes,r_mode)

    r_prompt_keymap = Dict{Any,Any}(
        '$' => function (s,args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, r_mode) do
                    LineEdit.state(s, r_mode).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s, '$')
            end
        end
    )

    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, r_prompt_keymap);
    nothing
end

function repl_inited(repl)
    mirepl = isdefined(repl,:mi) ? repl.mi : repl
    any(:prompt in fieldnames(typeof(m)) && m.prompt == "R> " for m in mirepl.interface.modes)
end

end # module


text:
using Statistics,Serialization,StatsPlots, DataFrames, LaTeXStrings

### baseline scenario in Germany with figure 1 a,b,c,d
### can only be plotted with green lines if empirical data is there
include("covid_plot_baseline_GER.jl")

### comparison plots w fig 2, supp fig 3
include("covid_plot_comparison_xi06.jl")
include("covid_plot_comparison_xi05.jl")

### and with dynamics for fig 3
include("covid_plot_dynamics_xi06.jl")
include("covid_plot_dynamics_xi05.jl")

### comparison plots including bailout porgram and shorttime program, supp fig 2
include("covid_plot_comparison_xi06_bailout.jl")
include("covid_plot_comparison_xi05_bailout.jl")

### no policy scenarios for supp fig 1
include("covid_plot_no_pol.jl")


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

# NDArray reshape (#272)
@deprecate reshape(arr::NDArray; shape=()) reshape(arr, shape)
@deprecate Reshape(arr::NDArray; shape=()) reshape(arr, shape)

# SymbolicNode reshape (#279)
@deprecate reshape(sym::SymbolicNode; shape=()) reshape(sym, shape)
@deprecate Reshape(sym::SymbolicNode; shape=()) reshape(sym, shape)

# srand (#282)
@deprecate srand!(seed_state::Int) srand(seed_state)

# v0.4
@deprecate sin(x::NDArray)    sin.(x)
@deprecate cos(x::NDArray)    cos.(x)
@deprecate tan(x::NDArray)    tan.(x)
@deprecate arcsin(x::NDArray) asin.(x)
@deprecate arccos(x::NDArray) acos.(x)
@deprecate arctan(x::NDArray) atan.(x)

@deprecate sinh(x::NDArray)    sinh.(x)
@deprecate cosh(x::NDArray)    cosh.(x)
@deprecate tanh(x::NDArray)    tanh.(x)
@deprecate arcsinh(x::NDArray) asinh.(x)
@deprecate arccosh(x::NDArray) acosh.(x)
@deprecate arctanh(x::NDArray) atanh.(x)

# @deprecate make `randn` exported accidentially
# so we make the depwarn manually
function randn(μ, σ, dims::NTuple{N,Int}, ctx::Context = cpu()) where N
  @warn("mx.randn(μ, σ, dims, ctx = cpu()) is deprecated, use " *
        "mx.randn(dims...; μ = μ, σ = σ, context = ctx) instead.")
  mx.randn(dims...; μ = μ, σ = σ, context = ctx)
end

function randn!(μ, σ, x::NDArray)
  @warn("mx.randn!(μ, σ, x::NDArray) is deprecated, use " *
        "mx.randn!(x; μ = μ, σ = σ) instead.")
  randn!(x; μ = μ, σ = σ)
end

function rand!(low::Real, high::Real, x::NDArray)
  @warn("rand!(low, high, x::NDArray) is deprecated, use " *
        "rand!(x, low = low, high = high) instead.")
  rand!(x, low = low, high = high)
end

function rand(low::Real, high::Real, dims::NTuple{N,Int}, context::Context = cpu()) where N
  @warn("rand!(low, high, dims, x::NDArray, context = cpu()) is deprecated, use " *
        "rand!(dims..., x; low = low, high = high, context = cpu()) instead.")
  rand(dims...; low = low, high = high, context = context)
end

@deprecate sigmoid(x::NDArray)                      sigmoid.(x)
@deprecate relu(x::NDArray)                         relu.(x)
@deprecate softmax(x::NDArray; axis = ndims(x))     softmax.(x, axis)
@deprecate log_softmax(x::NDArray; axis = ndims(x)) log_softmax.(x, axis)

@deprecate clip(x; a_min = 0, a_max = 0) clip(x, a_min, a_max)

function broadcast_plus(x::NDArray, y::NDArray)
  @warn("broadcast_plus(x, y) is deprecated, use x .+ y instead.")
  x .+ y
end

function broadcast_add(x::NDArray, y::NDArray)
  @warn("broadcast_add(x, y) is deprecated, use x .+ y instead.")
  x .+ y
end

function broadcast_sub(x::NDArray, y::NDArray)
  @warn("broadcast_sub(x, y) is deprecated, use x .- y instead.")
  x .- y
end

function broadcast_minus(x::NDArray, y::NDArray)
  @warn("broadcast_minus(x, y) is deprecated, use x .- y instead.")
  x .- y
end

function broadcast_mul(x::NDArray, y::NDArray)
  @warn("broadcast_mul(x, y) is deprecated, use x .* y instead.")
  x .* y
end

function broadcast_div(x::NDArray, y::NDArray)
  @warn("broadcast_div(x, y) is deprecated, use x ./ y instead.")
  x ./ y
end

function broadcast_mod(x::NDArray, y::NDArray)
  @warn("broadcast_mod(x, y) is deprecated, use x .% y instead.")
  x .% y
end

function broadcast_power(x::NDArray, y::NDArray)
  @warn("broadcast_power(x, y) is deprecated, use x.^y instead.")
  x.^y
end

function broadcast_equal(x::NDArray, y::NDArray)
  @warn("broadcast_equal(x, y) is deprecated, use x .== y instead.")
  x .== y
end

function broadcast_not_equal(x::NDArray, y::NDArray)
  @warn("broadcast_not_equal(x, y) is deprecated, use x .== y instead.")
  x .!= y
end

function broadcast_greater(x::NDArray, y::NDArray)
  @warn("broadcast_greater(x, y) is deprecated, use x .== y instead.")
  x .> y
end

function broadcast_greater_equal(x::NDArray, y::NDArray)
  @warn("broadcast_greater_equal(x, y) is deprecated, use x .== y instead.")
  x .>= y
end

function broadcast_lesser(x::NDArray, y::NDArray)
  @warn("broadcast_lesser(x, y) is deprecated, use x .== y instead.")
  x .< y
end

function broadcast_lesser_equal(x::NDArray, y::NDArray)
  @warn("broadcast_lesser_equal(x, y) is deprecated, use x .== y instead.")
  x .<= y
end

function broadcast_maximum(x::NDArray, y::NDArray)
  @warn("broadcast_maximum(x, y) is deprecated, use max.(x, y) instead.")
  max.(x, y)
end

function broadcast_minimum(x::NDArray, y::NDArray)
  @warn("broadcast_minimum(x, y) is deprecated, use min.(x, y) instead.")
  min.(x, y)
end

function broadcast_hypot(x::NDArray, y::NDArray)
  @warn("broadcast_hypot(x, y) is deprecated, use hypot.(x, y) instead.")
  hypot.(x, y)
end

# Introduced by https://github.com/apache/incubator-mxnet/pull/12845
import Base: sum, maximum, minimum, prod, cat
@deprecate sum(x::NDArray, dims) sum(x, dims = dims)
@deprecate maximum(x::NDArray, dims) maximum(x, dims = dims)
@deprecate minimum(x::NDArray, dims) minimum(x, dims = dims)
@deprecate prod(x::NDArray, dims) prod(x, dims = dims)
@deprecate cat(dims, As::NDArray{T}...) where T cat(As..., dims = dims)

import Statistics: mean
@deprecate mean(x::NDArray, dims) mean(x, dims = dims)

max_stars_repo_name: hsali/incubator-mxnet

max_issues_repo_name: hsali/incubator-mxnet

max_forks_repo_name: hsali/incubator-mxnet


text:
using DifferentialEquations
using Plots
plotlyjs()

function ex2_case_2(IC)

#constants
g=9.81        #acceleration due to gravity [ms^-2]
m=3000        #mass of capsule [Kg]
L=10          #height of A-frame [m]

#initial conditions
r₀,v₀,u₀,traj=IC #traj-trajectory of capsule till now-not needed for calculation
ζ₀=π/6 # angle subtended by the A-frame [rad]
β₀=dζ=90*π/(180*45)

#function to compute lagrangian
function lagrangian!(du,u,p,t)
    θ,ω,ϕ,ψ=u
    r=r₀
    β=dζ=β₀
    ζ=ζ₀+β*t
    du[1]=dθ=ω
    du[2]=dω=(ψ^2)*sin(θ)*cos(θ) + (g*sin(θ)/r) -
        (L*(β^2)*(cos(ζ)*cos(θ)*cos(ϕ) + sin(ζ)*sin(θ))/r)
    du[3]=dϕ=ψ
    du[4]=dψ=(L*(β^2)*cos(ζ)*sin(ϕ)/(r*sin(θ))) - (2*ω*ψ*cot(θ))
end

#time interval
tspan=(0.0,45.0)

#solving
prob=ODEProblem(lagrangian!,u₀,tspan)
sol=solve(prob,saveat=0.1)

#post-processing
U=sol[1:end,:]
(θ,ω,ϕ,ψ)=[U[x,:] for x in 1:size(U,1)]
t=sol.t
r=5 .* ones(size(t))
β=dζ=β₀
ζ=ζ₀.+β.*t
x=r.*sin.(θ).*cos.(ϕ)-L.*cos.(ζ)
y=r.*sin.(θ).*sin.(ϕ)
z=r.*cos.(θ)+L.*sin.(ζ)
uₜ=U[:,end]
rₜ=r[end]
vₜ=-v₀

# changing coordinate frames (no change)
trajectory=hcat(x,y,z)
FC=(rₜ,vₜ,uₜ,trajectory)

#plotting the solution
 display(plot(x,y,z,label="trajectory",
    xaxis="x (m)",yaxis="y (m)",zaxis="z (m)",
    grid=(:on,:black)))
savefig("./plots/trajectory-2.png")

 p1=plot(sol,vars=(0,1),label="θ",
    color="green", xaxis="t (s)", yaxis="θ (rad)")
 p2=plot(sol,vars=(0,2),label="ω",
    color="orange", xaxis="t (s)", yaxis="ω (rad/s)")
 p3=plot(sol,vars=(0,3),label="ϕ",
    color="purple", xaxis="t (s)", yaxis="ϕ (rad)")
 p4=plot(sol,vars=(0,4),label="ψ",
    color="blue", xaxis="t (s)", yaxis="ψ (rad/s)")
display(plot(p1,p2,p3,p4,layout=(2,2)))
savefig("./plots/plots-2.png")

#energy plot
KE=0.5*m*((L.*β).^2 .+ (r.*ω).^2 .+ (r.*ψ.*sin.(θ)).^2 .+
     2*L.*r.*β.*(ω.*sin.(ζ).*cos.(θ).*cos.(ϕ)
     .-ω.*cos.(ζ).*sin.(θ).-β.*sin.(θ).*sin.(ζ).*sin.(ϕ)))
PE=m*g.*z
TE=KE.+PE
p7=plot(t,KE,label="KE",color="black",
    xaxis="t (s)", yaxis="KE (J)")
p8=plot(t,PE,label="PE",color="green",
    xaxis="t (s)", yaxis="PE (J)")
p9=plot(t,TE,label="TE",color="red",
    xaxis="t (s)", yaxis="TE (J)")
display(plot!(p7,p8,p9))
savefig("./plots/energyplot-2.png")

return FC

end

max_stars_repo_path: Modeling the Dragon Module as a Pendulum/Scripts/ex2-case-2.jl


text:
abstract type LevyProcess{T} <: ContinuousTimeProcess{T} end

"""
    CompoundPoisson{T} <: LevyProcess{T}

Abstract type. For a compound Poisson process define `rjumpsize(P) -> T` and
 `arrival(P) -> Distribution`.
"""
abstract type CompoundPoisson{T} <: LevyProcess{T} end

"""
    ExpCounting(λ)

Counting process with arrival times `arrival(P) = Exponential(1/λ)` and unit jumps.
"""
struct ExpCounting <: CompoundPoisson{Int}
    λ::Float64
end

arrival(P::ExpCounting) = Exponential(inv(P.λ))

function rjumpsize(::ExpCounting)
    1
end

function sample(T, P::CompoundPoisson)
    dt = rand(arrival(P))
    t = 0.0
    y = 0.0
    tt = [t]
    yy = [y]

    while t + dt <= T
        t = t + dt
        y = y + rjumpsize(P)
        append!(tt, t)
        append!(yy, y)

        dt = rand(arrival(P))
    end
    SamplePath(tt, yy)
end

"""
    GammaProcess

A *GammaProcess* with jump rate `γ` and inverse jump size `λ` has increments `Gamma(t*γ, 1/λ)` and Levy measure

```math
ν(x)=γ x^{-1}\\exp(-λ x),
```

Here `Gamma(α,θ)` is the Gamma distribution in julia's parametrization with shape parameter `α` and scale `θ`.
"""
struct GammaProcess <: LevyProcess{Float64}
    γ::Float64
    λ::Float64
end

"""
    uniform_thinning!(X, P::GammaProcess, γᵒ)

Return a Gamma process `Y` with new intensity `γᵒ`, such that
`X-Y` has intensity `γ-γᵒ` and `Y` and `X-Y` are independent.
In the limit ``dt \\to \\infty`` the new Gamma process has each of is jump removed with
probability `γᵒ/γ`. Overwrites `X` with `Y`.
"""
function uniform_thinning!(X, P::GammaProcess, γᵒ)
    tt = X.tt
    yy = X.yy
    γ = P.γ
    if γᵒ > γ
        throw(ArgumentError("γᵒ > γ"))
    end
    y = yy[1]
    for i in 2:length(tt)
        dt = tt[i] - tt[i-1]
        y, yy[i] = yy[i], yy[i-1] + (yy[i] - y)*rand(Beta(dt*γᵒ,  dt*(γ-γᵒ)))
    end
    X
end
struct VarianceGammaProcess <: LevyProcess{Float64}
    θ::Float64
    σ::Float64
    ν::Float64
end

struct VarianceGamma
    θ::Float64
    σ::Float64
    t::Float64
    ν::Float64
end

"""
    GammaBridge(t, v, P)

A `GammaProcess` `P` conditional on htting `v` at time `t`.
"""
struct GammaBridge  <: ContinuousTimeProcess{Float64}
    t::Float64
    v::Float64
    P::GammaProcess
end

sample(tt::AbstractVector{Float64}, P::LevyProcess{T}, x1=zero(T)) where {T} =
    sample!(SamplePath{T}(collect(tt), zeros(T, length(tt))), P, x1)

function sample!(X, P::LevyProcess{T}, x1=zero(T)) where {T}
    tt = X.tt
    yy = X.yy

    yy[1] = x = x1
    for i in 2:length(tt)
        x = x + rand(increment(tt[i]-tt[i-1], P))
        yy[i] = x
    end
    X
end

increment(t, P::GammaProcess) = Gamma(t*P.γ, 1/P.λ)

lp(s, x, t, y, P::GammaProcess) = logpdf(increment(t-s, P), y-x)

levy(P::GammaProcess, x) = P.γ/x*exp(-P.λ*x)
llevyx(P, x) = log(levy(P, x)*x)

increment(t, P::VarianceGammaProcess) = VarianceGamma(P.θ, P.σ, t, P.ν)

function rand(P::VarianceGamma)
    Z = randn()
    G = rand(Gamma(P.t/P.ν, P.ν))
    P.θ*G + P.σ*sqrt(G)*Z
end

function sample(tt::AbstractVector{Float64}, P::GammaBridge, x1::Float64 = 0.0)
    tt = collect(tt)
    t = P.t
    r = searchsorted(tt, t)
    if isempty(r) # add t between n = last(r) and first(r)=n+1
       tt = Float64[tt[1:last(r)]; t; tt[first(r):end]] # t now at first(r)
    end
    X = sample(tt, P.P, zero(x1))
    dx = P.v - x1
    yy = X.yy
    yy[:] =  yy .* (dx/yy[first(r)]) .+ x1
    if isempty(r) # remove (t,x1)
        tt = [tt[1:last(r)]; tt[first(r)+1:end]]
        yy = [yy[1:last(r)]; yy[first(r)+1:end]]
    end
    SamplePath{Float64}(tt, yy)
end

function sample!(X, P::GammaBridge, x1::Float64 = 0.0)
    tt = X.tt
    t = P.t
    r = searchsorted(tt, t)
    if isempty(r)
        throw(BoundsError("$tt $t"))
    end
    sample!(X, P.P, zero(x1))
    dx = P.v - x1
    yy = X.yy
    yy[:] =  yy .* (dx/yy[first(r)]) .+ x1
    X
end

"""
    LocalGammaProcess
"""
struct LocalGammaProcess
    P::GammaProcess
    θ::Vector{Float64} # slope parameter of length N
    ρ::Vector{Float64} # intercept parameter of length N
    b::Vector{Float64} # grid points b1, ..., bN
    LocalGammaProcess(P, θ, ρ, b) = length(θ) != length(b) ?
        throw(ArgumentError("θ and b differ in length")) : new(P, θ, ρ, b)

end

"""
    θ(x, P::LocalGammaProcess)

Inverse jump size compared to gamma process with same alpha and beta.
"""
function θ(x, P::LocalGammaProcess)
    N = length(P.θ)
    N == 0 && return 0.0
    x <= P.b[1] && return 0.0
    x > P.b[N] && return P.θ[N] * x + P.ρ[N]

    #k = Int(div(x, P.b1))
    k = first(searchsorted(P.b, x)) - 1
    P.θ[k] * x + P.ρ[k]
end

"""
     nu(k, P)

(Bin-wise) integral of the Levy measure ``\\nu(B_k)`` (sic).
"""
function nu(k, P)
    if k == 0 && length(P.b) == 0
        P.P.γ*(-log(P.P.λ))
    elseif k == 0
        P.P.γ*(-log(P.P.λ) - expint1((P.P.λ)*P.b[1])) # up to certain constant
    elseif k == length(P.θ)
        @assert((P.P.λ + P.θ[k]) > 0.0)
        P.P.γ*exp(-P.ρ[k])*(expint1((P.P.λ + P.θ[k])*P.b[k])) # - 0 (upper limit infty)
    else
        P.P.γ*exp(-P.ρ[k])*(expint1((P.P.λ + P.θ[k])*P.b[k]) - expint1((P.P.λ + P.θ[k])*P.b[k+1]))
    end
end

"""
    compensator(kstart, P::LocalGammaProcess)

Compensator of LocalGammaProcess
For `kstart = 1`, this is ``\\sum_{k=1}^N \\nu(B_k)``,
for `kstart = 0`, this is ``\\sum_{k=0}^N \\nu(B_k) - C`` (where ``C`` is a constant).
"""
function compensator(kstart, P::LocalGammaProcess)
    s = 0.0
    for k in kstart:length(P.θ)
        s = s + nu(k, P)
    end
    s
end

"""
    compensator0(kstart, P::LocalGammaProcess)

Compensator of GammaProcess approximating the LocalGammaProcess.
For `kstart == 1` (only choice) this is ``\\nu_0([b_1,\\infty)``.
"""
function compensator0(kstart, P::LocalGammaProcess)
    if kstart == 1
        return length(P.b) > 0.0 ? P.P.γ * (expint1(P.P.λ * P.b[1])) : 0.0
    else
        throw(ArgumentError("k != 1"))
    end
end

"""
    llikelihood(X::SamplePath, Pº::LocalGammaProcess, P::LocalGammaProcess)

Log-likelihood `dPº/dP`. (Up to proportionality.)
"""
function llikelihood(X::SamplePath, Pº::LocalGammaProcess, P::LocalGammaProcess, c = 0.0)::Float64
    @assert(Pº.P.γ == P.P.γ)
    if Pº.P.λ == P.P.λ # case 1: same on the first bin
        ll = 0.
        for i in 2:length(X.tt)
            dx = X.yy[i] - X.yy[i-1] - c
            ll = ll - (θ(dx, Pº) - θ(dx, P)) # θ(x, P) ≈ θ_k dx + ρ_k
        end
        ll = ll - (X.tt[end]-X.tt[1])*(compensator(1, Pº)-compensator(1, P))
        return ll
    elseif Pº.θ === P.θ && Pº.ρ === P.ρ # case 2:
        ll = -(Pº.P.λ - P.P.λ) * (X.yy[end] - X.yy[1])
        return ll - (X.tt[end] - X.tt[1])*(compensator(0, Pº)-compensator(0, P))
    else
        ll = 0.0
        u = X.yy[end] - X.yy[1]
        @assert(c == 0)
        for i in 2:length(X.tt)
            dx = X.yy[i] - X.yy[i-1]
            if dx > P.b[1]
                ll = ll - ((Pº.P.λ - P.P.λ)*dx + θ(dx, Pº) - θ(dx, P))
                u -= dx
            end
        end
        ll += -(Pº.P.λ - P.P.λ) * u
        ll += -(X.tt[end] - X.tt[1])*(compensator(0, Pº)-compensator(0, P))
        return ll
    end
end

"""
    llikelihood(X::SamplePath, P::LocalGammaProcess)

Bridge log-likelihood with respect to reference measure `P.P`.
(Up to proportionality.)
"""
function llikelihood(X::SamplePath, P::LocalGammaProcess, c = 0.0)::Float64
    ll = 0.
    for i in 2:length(X.tt)
        dx = X.yy[i] - X.yy[i-1] - c
        ll = ll - θ(dx, P)
    end
    ll = ll - (X.tt[end] - X.tt[1])*(compensator(1, P) - compensator0(1, P))
    return ll
end

function llikelihood(X::SamplePath, P::Union{LevyProcess,Wiener})
    ll = 0.0
    for i in 1:length(X)-1
        dt = X.tt[i+1]-X.tt[i]
        dx = X.yy[i+1]-X.yy[i]
        ll += logpdf(increment(dt, P), dx) #Normal(0, sqrt(dt))
    end
    ll
end

"""
    posterior(Val{:λ}, P::GammaProcess, U::SamplePath, prior = (0.0, 0.0))

Marginal posterior distribution of parameter `λ`. Interpretation of conjugate
prior is "observed time, observed increment".
"""
function posterior(::Val{:λ}, P::GammaProcess, U::SamplePath, prior = (0.0, 0.0))
     Gamma((prior[1] + U.tt[end] - U.tt[1])*P.γ, 1/(prior[2] + U.yy[end] - U.yy[1]))
end
#=
Note:

Prior for the rate of a Gamma distribution (corresponding to `λ` in GammaProcess)
is Gamma with rate "prior + total increment", but Julia parametrizes `Gamma` in scale.
=#


text:
using MLStyle: @match

replace_variable(ex::Expr, list::Dict) = replace_variable(Val(ex.head), ex, list)
replace_variable(ex::LineNumberNode, list::Dict) = ex
replace_variable(::Val, ex, list) = copy(ex) # skip other patterns

function replace_variable(::Val{:->}, ex::Expr, list::Dict)
    Expr(:->, ex.args[1], replace_variable(ex.args[2], list))
end

function replace_variable(::Val{:block}, ex::Expr, list::Dict)
    out = Expr(:block)
    for each in ex.args
        push!(out.args, replace_variable(each, list))
    end
    return out
end

function replace_variable(ex::Symbol, list::Dict)
    ex in keys(list) && return list[ex]
    return ex
end

function replace_variable(::Val{:kw}, ex::Expr, list::Dict)
    if ex.args[1] in keys(list)
        return Expr(:kw, ex.args[1], list[ex.args[1]])
    else
        return ex
    end
end

function replace_variable(::Val{:parameters}, ex::Expr, list::Dict)
    out = Expr(:parameters)
    for each in ex.args
        push!(out.args, replace_variable(each, list))
    end
    return out
end

function replace_variable(::Val{:call}, ex::Expr, list::Dict)
    out_ex = Expr(:call)
    for (k, each) in enumerate(ex.args)
        # expand ...
        e = @match each begin
            Expr(:..., x) =>
                if x in keys(list)
                    val = list[x]
                    append!(out_ex.args, val)
                else
                    push!(out_ex.args, each)
                end

            Expr(:parameters, xs...) => push!(out_ex.args, replace_variable(Val(:parameters), each, list))
            _ => push!(out_ex.args, replace_variable(each, list))
        end
    end

    return out_ex
end

replace_variable(x::Any, ::Dict) = x


text:
export Dirichlet

"""
Description:

    Dirichlet factor node

    Multivariate:
    f(out, a) = Dir(out|a)
              = Γ(Σ_i a_i)/(Π_i Γ(a_i)) Π_i out_i^{a_i}
    where 'a' is a vector with every a_i > 0

    Matrix variate:
    f(out, a) = Π_k Dir(out|a_*k)
    where 'a' represents a left-stochastic matrix with every a_jk > 0

Interfaces:

    1. out
    2. a

Construction:

    Dirichlet(id=:some_id)
"""
mutable struct Dirichlet <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function Dirichlet(out, a; id=generateId(Dirichlet))
        @ensureVariables(out, a)
        self = new(id, Array{Interface}(undef, 2), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:a] = self.interfaces[2] = associate!(Interface(self), a)

        return self
    end
end

slug(::Type{Dirichlet}) = "Dir"

format(dist::ProbabilityDistribution{V, Dirichlet}) where V<:VariateType = "$(slug(Dirichlet))(a=$(format(dist.params[:a])))"

ProbabilityDistribution(::Type{MatrixVariate}, ::Type{Dirichlet}; a=ones(3,3)) = ProbabilityDistribution{MatrixVariate, Dirichlet}(Dict(:a=>a))
ProbabilityDistribution(::Type{Multivariate}, ::Type{Dirichlet}; a=ones(3)) = ProbabilityDistribution{Multivariate, Dirichlet}(Dict(:a=>a))
ProbabilityDistribution(::Type{Dirichlet}; a=ones(3)) = ProbabilityDistribution{Multivariate, Dirichlet}(Dict(:a=>a))

dims(dist::ProbabilityDistribution{Multivariate, Dirichlet}) = length(dist.params[:a])
dims(dist::ProbabilityDistribution{MatrixVariate, Dirichlet}) = size(dist.params[:a])

vague(::Type{Dirichlet}, dims::Int64) = ProbabilityDistribution(Multivariate, Dirichlet, a=ones(dims))
vague(::Type{Dirichlet}, dims::Tuple{Int64}) = ProbabilityDistribution(Multivariate, Dirichlet, a=ones(dims))
vague(::Type{Dirichlet}, dims::Tuple{Int64, Int64}) = ProbabilityDistribution(MatrixVariate, Dirichlet, a=ones(dims))

isProper(dist::ProbabilityDistribution{V, Dirichlet}) where V<:VariateType = all(dist.params[:a] .> 0.0)

unsafeMean(dist::ProbabilityDistribution{Multivariate, Dirichlet}) = dist.params[:a]./sum(dist.params[:a])
unsafeMean(dist::ProbabilityDistribution{MatrixVariate, Dirichlet}) = dist.params[:a]./sum(dist.params[:a],dims=1) # Normalize columns

function unsafeLogMean(dist::ProbabilityDistribution{Multivariate, Dirichlet})
    a = clamp.(dist.params[:a], tiny, huge)
    return digamma.(a) .- digamma.(sum(a))
end
function unsafeLogMean(dist::ProbabilityDistribution{MatrixVariate, Dirichlet})
    a = clamp.(dist.params[:a], tiny, huge)
    return digamma.(a) .- digamma.(sum(a,dims=1)) # Normalize columns
end

logPdf(dist::ProbabilityDistribution{Multivariate, Dirichlet}, x) = sum((dist.params[:a].-1).*log.(x)) - sum(loggamma.(dist.params[:a])) + loggamma(sum(dist.params[:a]))
logPdf(dist::ProbabilityDistribution{MatrixVariate, Dirichlet}, x) = sum(sum((dist.params[:a].-1).*log.(x),dims=1) - sum(loggamma.(dist.params[:a]), dims=1) + loggamma.(sum(dist.params[:a],dims=1)))

function unsafeVar(dist::ProbabilityDistribution{Multivariate, Dirichlet})
    a_sum = sum(dist.params[:a])
    return dist.params[:a].*(a_sum .- dist.params[:a])./(a_sum^2*(a_sum + 1.0))
end

function prod!( x::ProbabilityDistribution{V, Dirichlet},
                y::ProbabilityDistribution{V, Dirichlet},
                z::ProbabilityDistribution{V, Dirichlet}=ProbabilityDistribution(V, Dirichlet, a=ones(dims(x)))) where V<:VariateType

    z.params[:a] = x.params[:a] + y.params[:a] .- 1.0

    return z
end

@symmetrical function prod!(x::ProbabilityDistribution{Multivariate, Dirichlet},
                            y::ProbabilityDistribution{Multivariate, PointMass},
                            z::ProbabilityDistribution{Multivariate, PointMass}=ProbabilityDistribution(Multivariate, PointMass, m=[NaN]))

    a_y = clamp.(y.params[:m], 0.0, 1.0)
    z.params[:m] = deepcopy(a_y)

    return z
end

@symmetrical function prod!(x::ProbabilityDistribution{MatrixVariate, Dirichlet},
                            y::ProbabilityDistribution{MatrixVariate, PointMass},
                            z::ProbabilityDistribution{MatrixVariate, PointMass}=ProbabilityDistribution(MatrixVariate, PointMass, m=mat(NaN)))

    A_y = clamp.(y.params[:m], 0.0, 1.0)
    z.params[:m] = deepcopy(A_y)

    return z
end

function sample(dist::ProbabilityDistribution{Multivariate, Dirichlet})
    smpl = Vector{Float64}(undef, length(dist.params[:a]))
    for (i, alpha_i) in enumerate(dist.params[:a])
        smpl[i] = sample(ProbabilityDistribution(Univariate,Gamma,a=alpha_i,b=1.0))
    end
    smpl = smpl./sum(smpl)
    return smpl
end

# Entropy functional
function differentialEntropy(dist::ProbabilityDistribution{Multivariate, Dirichlet})
    a_sum = sum(dist.params[:a])

    -sum( (dist.params[:a] .- 1.0).*(digamma.(dist.params[:a]) .- digamma.(a_sum)) ) -
    labsgamma(a_sum) +
    sum(labsgamma.(dist.params[:a]))
end

function differentialEntropy(dist::ProbabilityDistribution{MatrixVariate, Dirichlet})
    H = 0.0
    for k = 1:dims(dist)[2] # For all columns
        a_sum = sum(dist.params[:a][:,k])

        H += -sum( (dist.params[:a][:,k] .- 1.0).*(digamma.(dist.params[:a][:,k]) .- digamma.(a_sum)) ) -
        labsgamma(a_sum) +
        sum(labsgamma.(dist.params[:a][:,k]))
    end

    return H
end

# Average energy functional
function averageEnergy(::Type{Dirichlet}, marg_out::ProbabilityDistribution{Multivariate}, marg_a::ProbabilityDistribution{Multivariate, PointMass})
    a_sum = sum(marg_a.params[:m])

    -labsgamma(a_sum) +
    sum(labsgamma.(marg_a.params[:m])) -
    sum( (marg_a.params[:m] .- 1.0).*unsafeLogMean(marg_out) )
end

function averageEnergy(::Type{Dirichlet}, marg_out::ProbabilityDistribution{Multivariate}, marg_a::ProbabilityDistribution{Multivariate, SampleList})
    samples, weights = marg_a.params[:s], marg_a.params[:w]
    S = length(weights) #number of samples
    log_gamma_of_sum, sum_of_log_gamma = 0.0, 0.0

    for s=1:S
        log_gamma_of_sum += weights[s]*loggamma(sum(samples[s]))
        sum_of_log_gamma += weights[s]*sum(loggamma.(samples[s]))
    end

    return -log_gamma_of_sum + sum_of_log_gamma -
            sum((unsafeMeanVector(marg_a) .- 1.0).*unsafeLogMean(marg_out))
end

function averageEnergy(::Type{Dirichlet}, marg_out::ProbabilityDistribution{MatrixVariate}, marg_a::ProbabilityDistribution{MatrixVariate, PointMass})
    (dims(marg_out) == dims(marg_a)) || error("Distribution dimensions must agree")

    log_mean_marg_out = unsafeLogMean(marg_out)

    H = 0.0
    for k = 1:dims(marg_out)[2] # For all columns
        a_sum = sum(marg_a.params[:m][:,k])

        H += -labsgamma(a_sum) +
        sum(labsgamma.(marg_a.params[:m][:,k])) -
        sum( (marg_a.params[:m][:,k] .- 1.0).*log_mean_marg_out[:,k] )
    end

    return H
end