text:
function _MKLfactorize(KK, Rhs::Vector,ifree::Vector; verbose = false)

    A  = KK[ifree,ifree]
    B  = Rhs[ifree]
    # Initialize the PARDISO internal data structures.
    ps = MKLPardisoSolver()

    if verbose
        set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
    end

    # First set the matrix type to handle general real symmetric matrices
    set_matrixtype!(ps, Pardiso.REAL_SYM_POSDEF)
    # Initialize the default settings with the current matrix type
    pardisoinit(ps)
    set_nprocs!(ps, Threads.nthreads())

    fix_iparm!(ps, :T)
    # Get the correct matrix to be sent into the pardiso function.
    # :N for normal matrix, :T for transpose, :C for conjugate
    A_pardiso = get_matrix(ps, A, :T)

    # Analyze the matrix and compute a symbolic factorization.
    set_phase!(ps, Pardiso.ANALYSIS)
    pardiso(ps, A_pardiso, B)

    # Compute the numeric factorization.
    set_phase!(ps, Pardiso.NUM_FACT)
    pardiso(ps, A_pardiso, B)

    return ps, A_pardiso

end

function _MKLsolve!(T::Vector, A_pardiso, ps::MKLPardisoSolver, Rhs::Vector{Float64}, ifree::Vector{Int64})

    B = Rhs[ifree]
    X = similar(B)

    # Compute the solutions X using the symbolic factorization.
    set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
    pardiso(ps, X, A_pardiso, B)

    @turbo for i ∈ 1:length(ifree)
        T[ifree[i]] =  X[i]
    end

end

function _MKLrelease!(ps::MKLPardisoSolver)
    # Free the PARDISO data structures.
    set_phase!(ps, Pardiso.RELEASE_ALL)
    pardiso(ps)
end

function _MKLpardiso!(T::Vector{Float64}, KK:: SparseMatrixCSC,Rhs::Vector{Float64},ifree::Vector{Int64})

    A = KK[ifree,ifree]
    B = Rhs[ifree]
    X = similar(B)

    # Initialize the PARDISO internal data structures.
    # ps = PardisoSolver()
    ps = MKLPardisoSolver()
    # if verbose
    #     set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
    # end

    # First set the matrix type to handle general real symmetric matrices
    set_matrixtype!(ps, Pardiso.REAL_SYM_POSDEF)
    # Initialize the default settings with the current matrix type
    pardisoinit(ps)
    set_nprocs!(ps, Threads.nthreads())

    # Remember that we pass in a CSC matrix to Pardiso, so need
    # to set the transpose iparm option.
    fix_iparm!(ps, :T)
    # Get the correct matrix to be sent into the pardiso function.
    # :N for normal matrix, :T for transpose, :C for conjugate
    A_pardiso = get_matrix(ps, A, :T)

    T[ifree] .= solve(ps, A_pardiso, B)

    # solve!(ps,X, A, B)

    # # Analyze the matrix and compute a symbolic factorization.
    # set_phase!(ps, Pardiso.ANALYSIS)
    # # set_perm!(ps, randperm(n))
    # pardiso(ps, A_pardiso, B)
    # # @printf("The factors have %d nonzero entries.\n", get_iparm(ps, 18))

    # # Compute the numeric factorization.
    # set_phase!(ps, Pardiso.NUM_FACT)
    # pardiso(ps, A_pardiso, B)
    # # @printf("The matrix has %d positive and %d negative eigenvalues.\n",
    # #         get_iparm(ps, 22), get_iparm(ps, 23))

    # # Compute the solutions X using the symbolic factorization.
    # set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
    # pardiso(ps, X, A_pardiso, B)
    # # @printf("PARDISO performed %d iterative refinement steps.\n", get_iparm(ps, 7))

    # # solve!(ps,X, A_pardiso, B)
    # # Free the PARDISO data structures.
    # set_phase!(ps, Pardiso.RELEASE_ALL)
    # pardiso(ps)

    # @tturbo for i ∈ 1:length(ifree)
    #     T[ifree[i]] =  X[i]
    # end

end # END PARDISO SOLVER N.1

# PARDISO SOLVER N.2 ================================================================
function _MKLpardiso(A, B::Vector{Float64})

    # Allocate solution vector
    X  = similar(B)

    # Initialize the PARDISO internal data structures.
    # ps = PardisoSolver()
    ps = MKLPardisoSolver()
    verbose = false
    if verbose
        set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
    end

    # set_nprocs!(ps, 4) # Sets the number of threads to use
    # First set the matrix type to handle general real symmetric matrices
    set_matrixtype!(ps, Pardiso.REAL_SYM_POSDEF)
    # Initialize the default settings with the current matrix type
    pardisoinit(ps)
    # Remember that we pass in a CSC matrix to Pardiso, so need
    # to set the transpose iparm option.
    fix_iparm!(ps, :T)
    # Get the correct matrix to be sent into the pardiso function.
    # :N for normal matrix, :T for transpose, :C for conjugate
    A_pardiso = get_matrix(ps, A, :T)
    solve!(ps, X, A_pardiso, B)
    # Free the PARDISO data structures.
    set_phase!(ps, Pardiso.RELEASE_ALL)
    pardiso(ps)

   return X

end # END PARDISO SOLVER N.2

# PARDISO 6.2 N.1 ========================================================================
function _pardiso!(T::Vector{Float64},KK:: SparseMatrixCSC,Rhs::Vector{Float64},ifree::Vector{Int64})

    A  = KK[ifree,ifree]
    B  = Rhs[ifree]
    X  = similar(B)

    # Initialize the PARDISO internal data structures.
    ps      = PardisoSolver()
    verbose = false
    if verbose
        set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
    end

    # First set the matrix type to handle general real symmetric matrices
    set_matrixtype!(ps, Pardiso.REAL_SYM_POSDEF)

    # Initialize the default settings with the current matrix type
    pardisoinit(ps)

    # Get the correct matrix to be sent into the pardiso function.
    # :N for normal matrix, :T for transpose, :C for conjugate
    A_pardiso = get_matrix(ps, A, :T)

    # Analyze the matrix and compute a symbolic factorization.
    set_phase!(ps, Pardiso.ANALYSIS)
    pardiso(ps, A_pardiso, B)

    # Compute the numeric factorization.
    set_phase!(ps, Pardiso.NUM_FACT)
    pardiso(ps, A_pardiso, B)

    # Compute the solutions X using the symbolic factorization.
    set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
    pardiso(ps, X, A_pardiso, B)

    # Free the PARDISO data structures.
    set_phase!(ps, Pardiso.RELEASE_ALL)
    pardiso(ps)

    # Refill global solution array
    for i ∈ 1:length(ifree)
        T[ifree[i]] =  X[i]
    end

end # END PARDISO 6.2 N.1

# PARDISO 6.2 N.2 ========================================================================
function _pardiso(A:: SparseMatrixCSC,B::Vector{Float64})

    # Initialize the PARDISO internal data structures.
    ps      = PardisoSolver()
    verbose = false
    if verbose
        set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
    end

    # First set the matrix type to handle general real symmetric matrices
    set_matrixtype!(ps, Pardiso.REAL_SYM_POSDEF)

    # Initialize the default settings with the current matrix type
    pardisoinit(ps)

    # Get the correct matrix to be sent into the pardiso function.
    # :N for normal matrix, :T for transpose, :C for conjugate
    A_pardiso = get_matrix(ps, A, :T)

    # Analyze the matrix and compute a symbolic factorization.
    set_phase!(ps, Pardiso.ANALYSIS)
    pardiso(ps, A_pardiso, B)

    # Compute the numeric factorization.
    set_phase!(ps, Pardiso.NUM_FACT)
    pardiso(ps, A_pardiso, B)

    # Compute the solutions X using the symbolic factorization.
    set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
    X = similar(B)
    pardiso(ps, X, A_pardiso, B)

    # Free the PARDISO data structures.
    set_phase!(ps, Pardiso.RELEASE_ALL)
    pardiso(ps)

    return X

end

function _Cholesky!(T::Vector{Float64},KK:: SparseMatrixCSC,Rhs::Vector{Float64},ifree::Vector{Int64})
    Rhs_free    = view(Rhs,ifree)
    A           = KK[ifree,ifree]
    A           = factorize(A)
    T[ifree]    = A.L'\(A.L\Rhs_free[A.p])
end # END CHOLESKY SOLVER

"""
Solve Symmetric Positive Definite system of eqs using via CHOLESKY factorization and allocate solution
in the un-constrained elements of the solution array containing allocated boundary conditions. Return
solution T and factorization F
"""
@inline function _CholeskyFactorizationSolve(T::Vector{Float64},KK,Rhs::Vector{Float64},ifree::Vector{Int64})
    K = KK[ifree,ifree]
    # F = cholesky((K')')
    # F = factorize((K')')
    F = lu(K)
    _CholeskyWithFactorization!(T,F,Rhs,ifree)
    return T,F
end
"""
Solve Symmetric Positive Definite system of eqs using via CHOLESKY factorization. Return
solution T and factorization F
"""
@inline function _CholeskyFactorizationSolve(KK, Rhs::Vector{Float64})
    A = KK[ifree,ifree]
    # F = cholesky((A')')
    F = factorize((A')')
    T = F\Rhs[ifree]
    return T, F
end

"""
Solve Symmetric Positive Definite system of eqs using A GIVEN factorization and allocate solution
in the un-constrained elements of the solution array containing allocated boundary conditions. In
place substitution of input solution array T
"""
@inline function _CholeskyWithFactorization!(T::Vector{Float64}, A, Rhs, ifree::Vector{Int64})
    T[ifree] .= A\Rhs[ifree]
end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

struct UUID
    value::UInt128
end
UUID(u::NTuple{2, UInt64}) = UUID((UInt128(u[1]) << 64) | UInt128(u[2]))
UUID(u::NTuple{4, UInt32}) = UUID((UInt128(u[1]) << 96) | (UInt128(u[2]) << 64) |
                                  (UInt128(u[3]) << 32) | UInt128(u[4]))

function convert(::Type{NTuple{2, UInt64}}, uuid::UUID)
    uuid = uuid.value
    hi = UInt64((uuid >> 64) & 0xffffffffffffffff)
    lo = UInt64(uuid & 0xffffffffffffffff)
    return (hi, lo)
end

function convert(::Type{NTuple{4, UInt32}}, uuid::UUID)
    uuid = uuid.value
    hh = UInt32((uuid >> 96) & 0xffffffff)
    hl = UInt32((uuid >> 64) & 0xffffffff)
    lh = UInt32((uuid >> 32) & 0xffffffff)
    ll = UInt32(uuid & 0xffffffff)
    return (hh, hl, lh, ll)
end

UInt128(u::UUID) = u.value

let groupings = [1:8; 10:13; 15:18; 20:23; 25:36]
    global UUID
    function UUID(s::AbstractString)
        s = lowercase(s)

        if !occursin(r"^[0-9a-f]{8}(?:-[0-9a-f]{4}){3}-[0-9a-f]{12}$", s)
            throw(ArgumentError("Malformed UUID string: $(repr(s))"))
        end

        u = UInt128(0)
        for i in groupings
            u <<= 4
            d = s[i] - '0'
            u |= 0xf & (d - 39*(d > 9))
        end
        return UUID(u)
    end
end

let groupings = [36:-1:25; 23:-1:20; 18:-1:15; 13:-1:10; 8:-1:1]
    global string
    function string(u::UUID)
        u = u.value
        a = Base.StringVector(36)
        for i in groupings
            a[i] = hex_chars[1 + u & 0xf]
            u >>= 4
        end
        a[24] = a[19] = a[14] = a[9] = '-'
        return String(a)
    end
end

print(io::IO, u::UUID) = print(io, string(u))
show(io::IO, u::UUID) = print(io, "UUID(\"", u, "\")")

ext: jl

max_stars_repo_name: ninjin/julia

max_issues_repo_name: amilsted/julia

max_forks_repo_name: amilsted/julia


text:
function cosx(x, n)
    an = 1
    sn = 0
    m = 2
    sqrX = x^2
    for i in 0:(n+1)
        sn += an
        an = - an * sqrX /(m-1)/m
        m += 2
    end

    return sn
end


text:
using Tmatrix
using Zygote

function objective_function(r_array_input, θ_array_input)
    # println(r_array)
    # println("I am <<objective_function>>, and the type of r_array is $(typeof(r_array))")
    # println("r_array = $r_array")

    # quadrupling_mesh
    r_θ_array = Tmatrix.quadruple_mesh_density(r_array_input, θ_array_input)
    r_array = r_θ_array[:,1]
    θ_array = r_θ_array[:,2]
    println(r_array)
    println(θ_array)
    ϕ_array = zeros(size(θ_array))

    T = Tmatrix.T_matrix_SeparateRealImag_arbitrary_mesh(
        n_max, wl_or_freq_input, input_unit, Eps_r_r_1, Eps_r_i_1, Mu_r_r_1, Mu_r_i_1, Eps_r_r_2, Eps_r_i_2, Mu_r_r_2, Mu_r_i_2,
        (collect(r_array)), θ_array, ϕ_array,
        rotationally_symmetric, symmetric_about_plane_perpendicular_z, BigFloat_precision
    )
    k1_complex = Tmatrix.get_WaveVector(
        wl_or_freq_input;
        input_unit=input_unit,
        Eps_r=Complex(Eps_r_r_1, Eps_r_i_1),
        Mu_r=Complex(Mu_r_r_1, Mu_r_i_1)
    )
    # return abs(T[1]) # if we need to minimize the first element of T-matrix

    # return sum(abs.(T)) # if we need to minimize the sum of all real and imag

    return Tmatrix.get_OrentationAv_scattering_CrossSection_from_Tmatrix(T, k1_complex) # if we need to minimize scattering cross section

    """
    return -1 * Tmatrix.get_OrentationAv_emissivity_from_Tmatrix(
        T,
        Complex(k1_r, k1_i),
        Tmatrix.calculate_surface_area_of_axisymmetric_particle(r_array, θ_array)
    ) # if we need to minimize emissivity
    """
end

function ∂objective_function(r_array, θ_array)
    # println("I am <<∂objective_function>>, and the type of r_array is $(typeof(r_array))")
    # println("r_array = $r_array")
    return Zygote.gradient(r_array -> objective_function(r_array, θ_array), (collect(r_array)))[1]
end

wl_or_freq_input = 1e-6; input_unit = "m"
n_max = 1
Eps_r_r_1 = 1.0; Eps_r_i_1 = 0.0
Eps_r_r_2 = 1.5; Eps_r_i_2 = 0.01
Mu_r_r_1 = 1.0; Mu_r_i_1 = 0.0
Mu_r_r_2 = 1.0; Mu_r_i_2 = 0.0

n_θ_points = 10
θ_array = collect(LinRange(1e-6, pi - 1e-6, n_θ_points))
r_array = 0.5 * wl_or_freq_input * ones(size(θ_array))
rotationally_symmetric = true
symmetric_about_plane_perpendicular_z = false
BigFloat_precision = nothing

learning_rate = 0.5e-18
global loss_here = -1e6
global ∂loss = 0
println("Starting optimization for T-matrix ...")
global n_iteration = 0
loss_array = []
loss_old = 9999999999999
while (n_iteration < 50) # (loss_here < 0)
    n_iteration += 1
    r_array = r_array .- learning_rate .* ∂loss

    loss_here = objective_function(r_array, θ_array)
    ∂loss = ∂objective_function(r_array, θ_array)
    append!(loss_array, loss_here)

    println()
    println("iteration #$n_iteration: loss_here = $loss_here, ∂loss = $∂loss")
    println("r_array = $r_array")

    xyz = vcat(Tmatrix.convert_coordinates_Sph2Cart.(r_array, θ_array, zeros(size(r_array)))...)
    p1 = plot(xyz_initial[:,1], xyz_initial[:,3], aspect_ratio=:equal, label="initial particle")
    p1 = plot!(xyz[:,1], xyz[:,3], aspect_ratio=:equal, label="particle after $n_iteration iterations", title="scattering cross section = $loss_here m")
    p2 = plot(1:length(loss_array), loss_array, xlabel="iteration #", ylabel="scattering cross section (m)")
    fig = plot(p1, p2, layout=(1, 2), size=(1200, 800))
    savefig(
        fig,
        "cache/iteration_particle_plots/maximizing_emissivity_07_12_2021__/particle_geom_iteration_$n_iteration.png"
    )
    # if (loss_here - loss_old) < 1e-6; break; end
    loss_old = loss_here
end


text:
using RecipesBase, Statistics, DimensionalData

@recipe f(::Type{Filterbank}, fb::Filterbank) = fb.data

@userplot Waterfall

@recipe function f(h::Waterfall)
    if h.args[1] isa Filterbank
        data = h.args[1].data
    elseif h.args[1] isa DimArray
        data = h.args[1]
    else
        @error "Argument must be a filterbank or DimArray"
    end

    # set up the subplots

    fc := :viridis
    legend := false
    link := :both
    framestyle := [:none :axes :none]
    grid := false
    layout := @layout [topav              _
                       heatmap{0.9w,0.9h} _]

    # Main Waterfall
    @series begin
        seriestype := :heatmap
        subplot := 2
        data
    end

    # these are common to both average plots
    linecolor := :black
    seriestype := :line

    # upper
    @series begin
        title := ""
        subplot := 1
        dropdims(mean(data,dims=Freq),dims=Freq)
    end
end


text:
module CompStat

using Plots
pyplot()

using LaTeXStrings
using Roots
using NLsolve
using Distributions
using Macros
using CustomTypes

export plot_cmpstat_bw, plot_cmpstat_colors, plot_transitions, @spec_funs, @dual_to_dual, @dual_to_dual_wages, fill_DualSS!

function max_zeros(z)
    if isempty(z)
        return 0.
    else
        return maximum(z)
    end
end

@def spec_funs begin
    Ez = exp(α^2/2)
    G(x) = cdf(LogNormal(0,α),x)
    g(x) = pdf(LogNormal(0,α),x)
    I(x) = exp(α^2/2)*cdf(Normal(0,1),((α^2-log(max(x,0)))/α))
    I0(x) = I(x) - x*(1-G(x))
end

@def classic_to_classic begin
    α₀ = (G(max(zp,zp_old))-G(zp_old))/(1-G(zp_old))
    β₀ = (G(max(zc, zc_old))-G(zc_old))/(1-G(zc_old))
    c₀ᶜ = (1-α₀)*nᶜ_old
    c₀⁰ = (1-β₀)*n⁰_old
    ω = -(s+λ)

    ρ₀ᶜ = (c₀ᶜ*μ + c₀⁰*(-ξ - ω) + nᶜ*(-μ - ξ))/(μ + ξ + ω)
    ρ₁ᶜ = c₀ᶜ - nᶜ - ρ₀ᶜ

    ρ₀⁰ = -(c₀ᶜ*μ + c₀⁰*(-ξ - ω) + n⁰*μ + n⁰*ξ - μ)/(μ + ξ + ω)
    ρ₁⁰ = c₀⁰ - n⁰ - ρ₀⁰

    np_c(t)  = nᶜ + ρ₀ᶜ*exp(ω*t) + ρ₁ᶜ*exp(-(μ+ξ)*t)
    np_0(t)  = n⁰ + ρ₀⁰*exp(ω*t) + ρ₁⁰*exp(-(μ+ξ)*t)
    nᵖ(t) = np_c(t) + np_0(t)
end

@def classic_to_dual begin

    α₁ = -δ/2 - μᵖ/2 - μᶠ/2 - ξ/2 + sqrt(δ^2 - 2*δ*μᵖ + 2*δ*μᶠ - 2*δ*ξ + μᵖ^2 + 2*μᵖ*μᶠ + 2*μᵖ*ξ + μᶠ^2 - 2*μᶠ*ξ + ξ^2)/2
    α₂ = -δ/2 - μᵖ/2 - μᶠ/2 - ξ/2 - sqrt(δ^2 - 2*δ*μᵖ + 2*δ*μᶠ - 2*δ*ξ + μᵖ^2 + 2*μᵖ*μᶠ + 2*μᵖ*ξ + μᶠ^2 - 2*μᶠ*ξ + ξ^2)/2
    α₀ = (G(max(zp, zp_old))-G(zp_old))/(1-G(zp_old))
    β₀ = (G(max(zc, zc_old))-G(zc_old))/(1-G(zc_old))
    c₀ᶜ = (1-α₀)*nᶜ_old
    c₀⁰ = (1-β₀)*n⁰_old
    ω = -(s+λ)

    ρ₀ᶜ = (c₀ᶜ*(δ*μᵖ + μᵖ*ω) + c₀⁰*(-δ*ξ - δ*ω - μᶠ*ξ - μᶠ*ω - ξ*ω - ω^2) - npᶜ*α₁*α₂ + nf_old*(μᵖ*ξ + μᵖ*ω) - μᵖ*ξ - μᵖ*ω)/(α₁*α₂ - α₁*ω - α₂*ω + ω^2)
    ρ₁ᶜ = (c₀ᶜ*(α₁ + δ + μᵖ + μᶠ) + c₀⁰*(-ξ - ω) + npᶜ*α₂ + ρ₀ᶜ*(α₂ - ω))/(α₁ - α₂)
    ρ₂ᶜ = c₀ᶜ - npᶜ - ρ₀ᶜ - ρ₁ᶜ

    ρ₀⁰ = -(c₀ᶜ*(δ*μᵖ + μᵖ*ω) + c₀⁰*(-δ*ξ - δ*ω - μᶠ*ξ - μᶠ*ω - ξ*ω - ω^2) + np⁰*α₁*α₂ + nf_old*(μᵖ*ξ + μᵖ*ω) - δ*μᵖ - μᵖ*ξ - μᵖ*ω)/(α₁*α₂ - α₁*ω - α₂*ω + ω^2)
    ρ₁⁰ = -(c₀ᶜ*μᵖ + c₀⁰*(-α₁ - δ - μᶠ - ξ - ω) - np⁰*α₂ + nf_old*μᵖ - μᵖ + ρ₀⁰*(-α₂ + ω))/(α₁ - α₂)
    ρ₂⁰ = c₀⁰ - np⁰ - ρ₀⁰ - ρ₁⁰

    ρ₁ᶠ = -(c₀ᶜ*μᶠ + c₀⁰*μᶠ - nf*α₂ + nf_old*(-α₁ - μᵖ - ξ) - μᶠ)/(α₁ - α₂)
    ρ₂ᶠ = -nf + nf_old - ρ₁ᶠ

    ρ₁ᵖ = ρ₁ᶜ + ρ₁⁰
    ρ₂ᵖ = ρ₂ᶜ + ρ₂⁰

    np_c(t)  = npᶜ + ρ₀ᶜ*exp(ω*t) + ρ₁ᶜ*exp(α₁*t) +  ρ₂ᶜ*exp(α₂*t)
    np_0(t)  = np⁰ + ρ₀⁰*exp(ω*t) + ρ₁⁰*exp(α₁*t) +  ρ₂⁰*exp(α₂*t)
    nᵖ(t) = np_c(t) + np_0(t)
    nᶠ(t) = nf + ρ₁ᶠ*exp(α₁*t) + ρ₂ᶠ*exp(α₂*t)
    u(t) = 1 - nᵖ(t) - nᶠ(t)
end

@def dual_to_dual begin

    α₁ = -δ/2 - μᵖ/2 - μᶠ/2 - ξ/2 + sqrt(δ^2 - 2*δ*μᵖ + 2*δ*μᶠ - 2*δ*ξ + μᵖ^2 + 2*μᵖ*μᶠ + 2*μᵖ*ξ + μᶠ^2 - 2*μᶠ*ξ + ξ^2)/2
    α₂ = -δ/2 - μᵖ/2 - μᶠ/2 - ξ/2 - sqrt(δ^2 - 2*δ*μᵖ + 2*δ*μᶠ - 2*δ*ξ + μᵖ^2 + 2*μᵖ*μᶠ + 2*μᵖ*ξ + μᶠ^2 - 2*μᶠ*ξ + ξ^2)/2
    α₀ = (G(max(zp,zp_old))-G(zp_old))/(1-G(zp_old))
    β₀ = (G(max(zc, zc_old, zs_old))-G(max(zc_old, zs_old)))/(1-G(max(zc_old, zs_old)))
    c₀ᶜ = (1-α₀)*npᶜ_old
    c₀⁰ = (1-β₀)*np⁰_old
    ω = -(s+λ)

    ρ₀ᶜ = (c₀ᶜ*(δ*μᵖ + μᵖ*ω) + c₀⁰*(-δ*ξ - δ*ω - μᶠ*ξ - μᶠ*ω - ξ*ω - ω^2) - npᶜ*α₁*α₂ + nf_old*(μᵖ*ξ + μᵖ*ω) - μᵖ*ξ - μᵖ*ω)/(α₁*α₂ - α₁*ω - α₂*ω + ω^2)
    ρ₁ᶜ = (c₀ᶜ*(α₁ + δ + μᵖ + μᶠ) + c₀⁰*(-ξ - ω) + npᶜ*α₂ + ρ₀ᶜ*(α₂ - ω))/(α₁ - α₂)
    ρ₂ᶜ = c₀ᶜ - npᶜ - ρ₀ᶜ - ρ₁ᶜ

    ρ₀⁰ = -(c₀ᶜ*(δ*μᵖ + μᵖ*ω) + c₀⁰*(-δ*ξ - δ*ω - μᶠ*ξ - μᶠ*ω - ξ*ω - ω^2) + np⁰*α₁*α₂ + nf_old*(μᵖ*ξ + μᵖ*ω) - δ*μᵖ - μᵖ*ξ - μᵖ*ω)/(α₁*α₂ - α₁*ω - α₂*ω + ω^2)
    ρ₁⁰ = -(c₀ᶜ*μᵖ + c₀⁰*(-α₁ - δ - μᶠ - ξ - ω) - np⁰*α₂ + nf_old*μᵖ - μᵖ + ρ₀⁰*(-α₂ + ω))/(α₁ - α₂)
    ρ₂⁰ = c₀⁰ - np⁰ - ρ₀⁰ - ρ₁⁰

    ρ₁ᶠ = -(c₀ᶜ*μᶠ + c₀⁰*μᶠ - nf*α₂ + nf_old*(-α₁ - μᵖ - ξ) - μᶠ)/(α₁ - α₂)
    ρ₂ᶠ = -nf + nf_old - ρ₁ᶠ

    ρ₁ᵖ = ρ₁ᶜ + ρ₁⁰
    ρ₂ᵖ = ρ₂ᶜ + ρ₂⁰

    np_c(t)  = npᶜ + ρ₀ᶜ*exp(ω*t) + ρ₁ᶜ*exp(α₁*t) +  ρ₂ᶜ*exp(α₂*t)
    np_0(t)  = np⁰ + ρ₀⁰*exp(ω*t) + ρ₁⁰*exp(α₁*t) +  ρ₂⁰*exp(α₂*t)
    nᵖ(t) = np_c(t) + np_0(t)
    nᶠ(t) = nf + ρ₁ᶠ*exp(α₁*t) + ρ₂ᶠ*exp(α₂*t)
    u(t) = 1 - nᵖ(t) - nᶠ(t)
end

@def dual_to_dual_wages begin

    Eᶜ = η*(I(zp)/(1-G(zp)) + (r+s+ϵ)*F - ϵ*F_bar + γ*θ) + (1-η)*b
    E₀ᶜ = 0.
    if zp > zp_old
        E₀ᶜ = η*((I(zp_old) - I(zp))/(G(zp) - G(zp_old)) + (r+s+ϵ)*F - ϵ*F_bar + γ*θ) + (1-η)*b
    end
    E⁰ = η*(I(zs)/(1-G(zs)) - λ*F + γ*θ) + (1-η)*b
    E₀⁰ = 0.
    if zc > zs_old
        E₀⁰ = η*((I(zs_old) - I(zc))/(G(zc)-G(zs_old)) + γ*θ) + (1-η)*b
    end

    Wpᶜ = (Eᶜ*npᶜ*ξ + Eᶜ*npᶜ*ω + Eᶜ*np⁰*ξ + Eᶜ*np⁰*ω)/ω
    ν₀ᶜ = -Eᶜ*ξ*ρ₀ᶜ - Eᶜ*ξ*ρ₀⁰ - Eᶜ*ρ₀ᶜ*ω - Eᶜ*ρ₀⁰*ω
    ν₁ᶜ = -Eᶜ*npᶜ - Eᶜ*np⁰ - Eᶜ*ξ*(npᶜ + np⁰)/ω - Eᶜ*ρ₁ᶜ - Eᶜ*ρ₁⁰ - Eᶜ*ρ₂ᶜ - Eᶜ*ρ₂⁰ + Eᶜ*(α₁ + ξ)*(ρ₁ᶜ + ρ₁⁰)/(α₁ - ω) + Eᶜ*(α₂ + ξ)*(ρ₂ᶜ + ρ₂⁰)/(α₂ - ω) - E₀ᶜ*npᶜ_old*α₀ + Wpᶜ_old
    ν₂ᶜ = (-Eᶜ*ξ*ρ₁ᶜ - Eᶜ*ξ*ρ₁⁰ - Eᶜ*ρ₁ᶜ*ω - Eᶜ*ρ₁⁰*ω)/(α₁ - ω)
    ν₃ᶜ = (-Eᶜ*ξ*ρ₂ᶜ - Eᶜ*ξ*ρ₂⁰ - Eᶜ*ρ₂ᶜ*ω - Eᶜ*ρ₂⁰*ω)/(α₂ - ω)

    Wp⁰ = (E⁰*npᶜ*μᵖ + E⁰*nf*μᵖ + E⁰*np⁰*μᵖ - E⁰*μᵖ)/ω
    ν₀⁰ = -E⁰*μᵖ*ρ₀ᶜ - E⁰*μᵖ*ρ₀⁰
    ν₁⁰ = E⁰*μᵖ*(ρ₂ᶜ + ρ₂ᶠ + ρ₂⁰)/(α₂ - ω) + E⁰*μᵖ*(ρ₁ᶜ + ρ₁ᶠ + ρ₁⁰)/(α₁ - ω) - E⁰*μᵖ*(npᶜ + nf + np⁰ - 1)/ω - E₀⁰*np⁰_old*β₀ + Wp⁰_old
    ν₂⁰ = (-E⁰*μᵖ*ρ₁ᶜ - E⁰*μᵖ*ρ₁ᶠ - E⁰*μᵖ*ρ₁⁰)/(α₁ - ω)
    ν₃⁰ = (-E⁰*μᵖ*ρ₂ᶜ - E⁰*μᵖ*ρ₂ᶠ - E⁰*μᵖ*ρ₂⁰)/(α₂ - ω)

    Wp_c(t) = Wpᶜ + ν₀ᶜ*t*exp(ω*t) + ν₁ᶜ*exp(ω*t) + ν₂ᶜ*exp(α₁*t) + ν₃ᶜ*exp(α₂*t)
    Wp_0(t) = Wp⁰ + ν₀⁰*t*exp(ω*t) + ν₁⁰*exp(ω*t) + ν₂⁰*exp(α₁*t) + ν₃⁰*exp(α₂*t)
    Wᵖ(t) = Wp_c(t) + Wp_0(t)

end

@def dual_to_classic begin
    α₀ = (G(max(zp,zp_old))-G(zp_old))/(1-G(zp_old))
    β₀ = (G(max(zc, zc_old))-G(zc_old))/(1-G(zc_old))
    c₀ᶜ = (1-α₀)*npᶜ_old
    c₀⁰ = (1-β₀)*np⁰_old
    ω = -(s+λ)

    ρ₀ᶜ = (c₀ᶜ*(δ*μ + μ*ω) + c₀⁰*(-δ*ξ - δ*ω - ξ*ω - ω^2) + nᶜ*(-δ*μ - δ*ξ) + nf_old*(μ*ξ + μ*ω) - μ*ξ - μ*ω)/(δ*μ + δ*ξ + δ*ω + μ*ω + ξ*ω + ω^2)
    ρ₁ᶜ = (c₀ᶜ*(δ - ξ) + c₀⁰*(-ξ - ω) - nᶜ*δ + ρ₀ᶜ*(-δ - ω))/(δ - μ - ξ)
    ρ₂ᶜ = c₀ᶜ - nᶜ - ρ₀ᶜ - ρ₁ᶜ

    ρ₀⁰ = -(c₀ᶜ*(δ*μ + μ*ω) + c₀⁰*(-δ*ξ - δ*ω - ξ*ω - ω^2) + n⁰*δ*μ + n⁰*δ*ξ + nf_old*(μ*ξ + μ*ω) - δ*μ - μ*ξ - μ*ω)/(δ*μ + δ*ξ + δ*ω + μ*ω + ξ*ω + ω^2)
    ρ₁⁰ = -(c₀ᶜ*μ + c₀⁰*(-δ + μ - ω) + n⁰*δ + nf_old*μ - μ + ρ₀⁰*(δ + ω))/(δ - μ - ξ)
    ρ₂⁰ = c₀⁰ - n⁰ - ρ₀⁰ - ρ₁⁰

    ρ₁ᵖ = ρ₁ᶜ + ρ₁⁰
    ρ₂ᵖ = ρ₂ᶜ + ρ₂⁰

    np_c(t)  = nᶜ + ρ₀ᶜ*exp(ω*t) + ρ₁ᶜ*exp(-(μ + ξ)*t) +  ρ₂ᶜ*exp(-δ*t)
    np_0(t)  = n⁰ + ρ₀⁰*exp(ω*t) + ρ₁⁰*exp(-(μ + ξ)*t) +  ρ₂⁰*exp(-δ*t)
    nᵖ(t) = np_c(t) + np_0(t)
    nᶠ(t) = nf_old*exp(-δ*t)
    u(t) = 1 - nᵖ(t) - nᶠ(t)

end

"""
    compute_dual!(ss, sol, param)

Fills the `DualSS` structure with the dual-equilibrium thresholds and employments given basic equilibrium values `sol` and parameters `param`
"""
function compute_dual!(ss, sol, param; index=nothing)

    Macros.@load_struct param η r s λ α m b γ σ F δ
    @spec_funs

    zp = sol[1]
    θ = sol[2]

    zc = zp+(r+s+λ)*F
    q = m*θ^(-σ)
    rU = b + η*γ*θ/(1-η)
    zf = rU - λ*(Ez-rU)/(r+δ)
    zs = (zp/(r+s+λ) + F - zf/(r+δ+λ))/(1/(r+s+λ) - 1/(r+δ+λ))
    ξ = s+λ*G(zp)
    p = m*θ^(1-σ)
    μᵖ = p*(1-G(max(zs,zc)))
    μᶠ = p*(G(max(zs,zf))-G(zf))
    npᶜ = λ*(1-G(zp))*μᵖ*δ/((s+λ)*(μᵖ*δ + ξ*δ + μᶠ*ξ))
    np⁰ = ξ*μᵖ*δ/((s+λ)*(μᵖ*δ + ξ*δ + μᶠ*ξ))
    np = npᶜ + np⁰
    nf = μᶠ*ξ/(μᵖ*δ + ξ*δ + μᶠ*ξ)
    u = 1-np-nf
    v = u*θ

    z_bar_p = (λ*np*I(zp)+p*u*I(max(zs,zc)))/((s+λ)*np)
    z_bar_f = (λ*nf*Ez+p*u*(I(zf)-I(max(zs,zf))))/((δ+λ)*nf)
    z_bar = (np*z_bar_p+nf*z_bar_f)/(np+nf)

    W = z_bar_p*np + z_bar_f*nf + b*u - γ*v - λ*G(zp)*np*F
    W /= r
    if isnothing(index)
        @save_struct ss zp zc zs zf npᶜ np⁰ np nf θ rU ξ μᵖ μᶠ z_bar_p z_bar_f z_bar W u
    else
        @save_struct_index ss index zp zc zs zf npᶜ np⁰ np nf θ rU ξ μᵖ μᶠ z_bar_p z_bar_f z_bar W u
    end

end

"""
    compute_classic!(ss, sol, param)

Fills the `ClassicSS` structure with the classic-equilibrium thresholds and employments given basic equilibrium values `sol` and parameters `param`
"""
function compute_classic!(ss, sol, param; index=nothing)
    Macros.@load_struct param η r s λ α m b γ σ F δ
    @spec_funs

    zp = sol[1]
    θ = sol[2]

    q = m*θ^(-σ)
    rU = b+η*γ*θ/(1-η)
    zc = zp+(r+s+λ)*F
    ξ = s+λ*G(zp)
    p = m*θ^(1-σ)
    μ = p*(1-G(zc))
    nᶜ = λ*(1-G(zp))*μ/((ξ+μ)*(s+λ))
    n⁰ = ξ*μ/((ξ+μ)*(s+λ))
    n = μ/(ξ+μ)
    u = 1-n
    v = u*θ

    z_bar = (λ*n*I(zp)+p*u*I(zc))/((s+λ)*n)
    W = z_bar*n + b*u - γ*v - λ*G(zp)*n*F
    W /= r
    if isnothing(index)
        @save_struct ss zp zc n nᶜ n⁰ θ μ z_bar W ξ rU u
    else
        @save_struct_index ss index zp zc n nᶜ n⁰ θ μ z_bar W ξ rU u
    end

end

"""
    EqClassic!(R,x,param)

Replaces `R` with the values of the fundamental equation block of the dual equilibrium given fundamental variables `x` and parameters `param`.
"""
function EqClassic!(R, x, param)
    Macros.@load_struct param η r s λ α m b γ σ F δ
    @spec_funs

    zp = x[1]
    θ = max(x[2],0.)

    q = m*θ^(-σ)
    rU = b + η*γ*θ/(1-η)
    zc = zp + (r+s+λ)*F

    R[1] = zp+(r+s)*F+λ*I0(zp)/(r+s+λ)-rU
    R[2] = γ - q*(1-η)*I0(zc)/(r+s+λ)

end

"""
    EqDual!(R,x,param)

Replaces `R` with the values of the fundamental equation block of the dual equilibrium given fundamental variables `x` and parameters `param`. Used to compute the steady state in the dual model.
"""
function EqDual!(R, x, param)
    Macros.@load_struct param η r s λ α m b γ σ F δ
    @spec_funs

    zp = x[1]
    θ = max(x[2],0.)

    zc = zp+(r+s+λ)*F
    q = m*θ^(-σ)
    rU = b + η*γ*θ/(1-η)
    zf = rU - λ*(Ez-rU)/(r+δ)
    zs = (zp/(r+s+λ) + F - zf/(r+δ+λ))/(1/(r+s+λ) - 1/(r+δ+λ))

    R[1] = zp+(r+s)*F+λ*I0(zp)/(r+s+λ)-rU
    R[2] = γ - q*(1-η)*( I0(max(zs,zc))/(r+s+λ) + (I0(zf)-I0(max(zs,zf)))/(r+δ+λ) )

end

"""
    solve_F1!(R, x, param)

Replaces `R` with the values of the fundamental equation block of the dual equilibrium under the constraint ``z^f = z^* = z^c``. Used to compute the dual steady state under the constraint ``z^f = z^* = z^c``.
"""
function solve_F1!(R, x, param)
    Macros.@load_struct param η r s λ α m b γ σ δ
    @spec_funs

    F1 = x[1]
    θ = max(x[2],0.)

    q = m*θ^(-σ)
    rU = b + η*γ*θ/(1-η)
    zc = rU - λ*(Ez-rU)/(r+δ)
    zp = zc - (r+s+λ)*F1

    R[1] = zp+(r+s)*F1+λ*I0(zp)/(r+s+λ)-rU
    R[2] = γ - q*(1-η)*I0(zc)/(r+s+λ)
end

"""
    solve_F2!(R, x, param)

Replaces `R` with the values of the fundamental equation block of the dual equilibrium under the constraint ``z^p = 0``. Used to compute the dual steady state under the constraint ``z^p = 0``.
"""
function solve_F2!(R, x, param)
    Macros.@load_struct param η r s λ α m b γ σ δ
    @spec_funs

    F2 = x[1]
    θ = max(x[2],0.)

    q = m*θ^(-σ)
    rU = b + η*γ*θ/(1-η)
    zf = rU - λ*(Ez-rU)/(r+δ)
    zs = ( F2 - zf/(r+δ+λ))/(1/(r+s+λ) - 1/(r+δ+λ))

    R[1] = (r+s)*F2+λ*Ez/(r+s+λ)-rU
    R[2] = γ - q*(1-η)*( I0(zs)/(r+s+λ) + (I0(zf)-I0(zs))/(r+δ+λ) )
end

"""
    fill_AdjDual!(adj_dual, F, ss_old, param, x0)

Fills the `index` index of the fields of the `comp_stat` structure with labor market values when firing costs equal `F`, `param` contains the other parameters, `ss` contains the initial steady state and `x0` is the initial guess for solving the system `EqDual` derives.
"""
function fill_AdjDual!(adj_dual, F, ss_old, param, x0; index=nothing)

    Macros.@load_struct param η r s λ α m b γ σ δ
    @spec_funs

    param.F = F
    ss = DualSS()
    fill_DualSS!(ss, param, x0)

    @load_struct_cat old ss_old npᶜ np⁰ np nf zp zc zs
    @load_struct ss zp zc zs zf npᶜ np⁰ np nf θ rU ξ μᵖ μᶠ z_bar_p z_bar_f z_bar W u

    @dual_to_dual

    τ_p = max_zeros(fzeros(t->abs((nᵖ(t)-np)/np_old)-0.01,0,500))
    τ_f = max_zeros(fzeros(t->abs((nᶠ(t)-nf)/nf_old)-0.01,0,500))

    if isnothing(index)
        @save_struct adj_dual τ_p τ_f
    else
        @save_struct_index adj_dual index τ_p τ_f
    end

end

"""
    fill_AdjClassic!(adj_classic, F, ss_old, param, x0)

Fills the fields of the `comp_stat` structure with labor market values when firing costs equal `F`, `param` contains the other parameters, `ss` contains the initial steady state and `x0` is the initial guess for solving the system `EqDual` derives.
"""
function fill_AdjClassic!(adj_classic, F, ss_old, param, x0; index=nothing)

    Macros.@load_struct param η r s λ α m b γ σ δ
    @spec_funs

    param.F = F
    ss = ClassicSS()
    fill_ClassicSS!(ss, param, x0)

    @load_struct_cat old ss_old nᶜ n⁰ n zp zc
    @load_struct ss zp zc n nᶜ n⁰ θ μ z_bar W ξ rU

    @classic_to_classic

    τ = max_zeros(fzeros(t->abs((nᵖ(t)-n)/n_old)-0.01,0,500))

    if isnothing(index)
        @save_struct adj_classic τ
    else
        @save_struct_index adj_classic index τ
    end

end

"""
    fill_ClassicSS!(ss, param, x0)

Fills the `ss` structure with the steady-state values of the classic model when `param` contains parameters and `x0` is the initial value to solve the system `EqClassic` derives.
"""
function fill_ClassicSS!(ss, param, x0; index=nothing)
    R = zeros(2)
    sol = nlsolve((R,x)->EqClassic!(R,x,param), x0, show_trace=false)
    compute_classic!(ss, sol.zero, param; index=index)
end

"""
    fill_DualSS!(ss, param, x0)

Fills the `ss` structure with the steady-state values of the dual model when `param` contains parameters and `x0` is the initial value to solve the system `EqDual` derives.
"""
function fill_DualSS!(ss, param, x0; index=nothing)
    R = zeros(2)
    sol = nlsolve((R,x)->EqDual!(R,x,param), x0, show_trace=false)
    compute_dual!(ss, sol.zero, param; index=index)
end

"""
    plot_cmpstat_bw(param, f; x0 = ..., x1 = ...)

Plots in black and white various comparative statics using parameters `param` and saving the figures with the suffix `f`. `x0` and `x1` contain starting values to solve the involved systems of equations.
"""
function plot_cmpstat_bw(param, f; x0 = [0.9, 1.3], x1 = [1.33, 1.3])

    @load_struct param F
    classic = ClassicSS()
    fill_ClassicSS!(classic, param, x0)
    dual = DualSS()
    fill_DualSS!(dual, param, x0)

    sol = nlsolve((R,x)->solve_F1!(R,x,param),x1,show_trace=true)
    F1 = sol.zero[1]
    sol = nlsolve((R,x)->solve_F2!(R,x,param),x1,show_trace=true)
    F2 = sol.zero[1]

    N_F = 100

    dual_new = DualSS()
    classic_new = ClassicSS()
    adj_dual = AdjDual()
    adj_classic = AdjClassic()
    init!(dual_new, N_F)
    init!(classic_new, N_F)
    init!(adj_dual, N_F)
    init!(adj_classic, N_F)

    grid_F = range(0.8*F1, stop=1.1*F, length=N_F)
    for i in 1:N_F
        param.F = grid_F[i]
        fill_ClassicSS!(classic_new, param, x0; index=i)
        fill_DualSS!(dual_new, param, x0; index=i)
        fill_AdjClassic!(adj_classic, grid_F[i], classic, param, x0; index=i)
        fill_AdjDual!(adj_dual, grid_F[i], dual, param, x0; index=i)
    end

    cd("figures")

    # Transition probabilities
    @load_struct dual_new ξ μᵖ μᶠ np nf u z_bar_p z_bar_f z_bar W
    @load_struct adj_dual τ_p τ_f
    @load_struct_cat c classic_new ξ μ n u z_bar W
    @load_struct_cat c adj_classic τ
    p1 = plot(grid_F, ξ,lab=L"$\xi$ - Dual",color=:black,line=(:solid, 1))
    plot!(grid_F, μᵖ,lab=L"$\mu^p$ - Dual",color=:black,line=(:dash, 2))
    plot!(grid_F, μᶠ,lab=L"$\mu^f$ - Dual",color=:black,line=(:dot, 1))
    plot!(grid_F, ξ_c,lab=L"$\xi$ - Classic",color=:black,line=(:dashdot, 1))
    plot!(grid_F, μ_c,lab=L"$\mu$ - Classic",color=:black,line=(:dash, 1), legend=:right)
    plot!([F1,F], seriestype=:vline, lab="",color=:black,line=(:dash, 1))
    xlabel!("Firing costs")
    ylabel!("Probabilities")
    xticks!([F1,F],[L"$\widehat{F}$",L"$F$"])
    xgrid!(:off)
    ygrid!(:off)
    savefig("transitions_"*f*".pdf")

    # Employment

    p2 = plot(grid_F, np,lab=L"$n^p$ - Dual",color=:black,line=(:dash, 2))
    plot!(grid_F, nf,lab=L"$n^f$ - Dual",color=:black,line=(:dot, 2))
    plot!(grid_F, u,lab=L"$u$ - Dual",color=:black,line=(:solid, 2))
    plot!(grid_F, n_c,lab=L"$n^p$ - Classic",color=:black,line=(:dash, 1))
    plot!(grid_F, u_c,lab=L"$u$ - Classic",color=:black,line=(:dashdot, 1), legend=:right)
    plot!([F1,F], seriestype=:vline, lab="",color=:black,line=(:dash, 1))
    xlabel!("Firing costs")
    ylabel!("Employment")
    xticks!([F1,F],[L"$\widehat{F}$",L"$F$"])
    xgrid!(:off)
    ygrid!(:off)
    savefig("employment_"*f*".pdf")

    #Probabilities and Employment
    plot(p1,p2,layout=(1,2))
    savefig("prob_emp_"*f*".pdf")

    # Productivity
    p4 = plot(grid_F, z_bar_p,lab=L"$\overline{z}^p$ - Dual",color=:black,line=(:dash, 2))
    plot!(grid_F, z_bar_f,lab=L"$\overline{z}^f$ - Dual",color=:black,line=(:dot, 2))
    plot!(grid_F, z_bar,lab=L"$\overline{z}$ - Dual",color=:black,line=(:solid, 2))
    plot!(grid_F, z_bar_c,lab=L"$\overline{z}$ - Classic",color=:black,line=(:dashdot, 1))
    plot!([F1,F], seriestype=:vline, lab="",color=:black,line=(:dash, 1))
    xlabel!("Firing costs")
    ylabel!("Productivity")
    xticks!([F1,F],[L"$\widehat{F}$",L"$F$"])
    xgrid!(:off)
    ygrid!(:off)
    savefig("productivity_"*f*".pdf")

    # Welfare
    p3 = plot(grid_F, W,lab=L"Dual",color=:black,line=(:solid, 2))
    plot!(grid_F, W_c,lab=L"Classic",color=:black,line=(:dash, 2))
    plot!([F1,F], seriestype=:vline, lab="",color=:black,line=(:dash, 1))
    xlabel!("Firing costs")
    ylabel!("Welfare")
    xticks!([F1,F],[L"$\widehat{F}$",L"$F$"])
    xgrid!(:off)
    ygrid!(:off)
    savefig("welfare_"*f*".pdf")

    #Welfare and Productivity
    plot(p3,p4,layout=(1,2))
    savefig("wel_prod_"*f*".pdf")

    # Adjustment
    plot(grid_F, τ_p,lab=L"$n^p$ - Dual",color=:black,line=(:solid, 2))
    plot!(grid_F, τ_f,lab=L"$n^f$ - Dual",color=:black,line=(:dot, 2))
    plot!(grid_F, τ_c,lab=L"$n$ - Classic",color=:black,line=(:dash, 2))
    plot!([F1,F], seriestype=:vline, lab="",color=:black,line=(:dash, 1))
    xlabel!("Firing costs")
    ylabel!(L"$\tau_{99\%}$")
    xticks!([F1,F],[L"$\widehat{F}$",L"$F$"])
    xgrid!(:off)
    ygrid!(:off)
    savefig("adjustment_"*f*".pdf")

    cd("..")

end

"""
    plot_cmpstat_colors(param, f; x0 = ..., x1 = ...)

Plots in colors various comparative statics using parameters `param` and saving the figures with the suffix `f`. `x0` and `x1` contain starting values to solve the involved systems of equations.
"""
function plot_cmpstat_colors(param, f; x0 = [0.9, 1.3], x1 = [1.33, 1.3])

    @load_struct param F
    classic = ClassicSS()
    fill_ClassicSS!(classic, param, x0)
    dual = DualSS()
    fill_DualSS!(dual, param, x0)

    sol = nlsolve((R,x)->solve_F1!(R,x,param),x1,show_trace=false)
    F1 = sol.zero[1]
    sol = nlsolve((R,x)->solve_F2!(R,x,param),x1,show_trace=false)
    F2 = sol.zero[1]

    N_F = 100

    dual_new = DualSS()
    classic_new = ClassicSS()
    adj_dual = AdjDual()
    adj_classic = AdjClassic()
    init!(dual_new, N_F)
    init!(classic_new, N_F)
    init!(adj_dual, N_F)
    init!(adj_classic, N_F)

    grid_F = range(0.8*F1, stop=1.1*F, length=N_F)
    for i in 1:N_F
        param.F = grid_F[i]
        fill_ClassicSS!(classic_new, param, x0; index=i)
        fill_DualSS!(dual_new, param, x0; index=i)
        fill_AdjClassic!(adj_classic, grid_F[i], classic, param, x0; index=i)
        fill_AdjDual!(adj_dual, grid_F[i], dual, param, x0; index=i)
    end

    cd("figures")

    # Transition probabilities
    @load_struct dual_new ξ μᵖ μᶠ np nf u z_bar_p z_bar_f z_bar W
    @load_struct adj_dual τ_p τ_f
    @load_struct_cat c classic_new ξ μ n u z_bar W
    @load_struct_cat c adj_classic τ
    p1 = plot(grid_F, ξ,lab="Open-ended JD",color=:red, line=(:dash, 1))
    plot!(grid_F, μᵖ,lab="Open-ended JC",color=:red, line=(:solid, 1))
    plot!(grid_F, μᶠ,lab="Fixed-term JC",color=:blue, line=(:solid, 1))
    plot!([F1,F], seriestype=:vline, lab="",color=:black, line=(:dash, 1))
    xlabel!("Firing costs")
    ylabel!("Probabilities")
    xticks!([F1,F],[L"$\widehat{F}$",L"$F$"])
    xgrid!(:off)
    ygrid!(:off)

    savefig("transitions_color_"*f*".pdf")

    # Employment
    p2 = plot(grid_F, nf, lab="Fixed-term",color=:blue,line=(:solid, 1))
    plot!([F1,F], seriestype=:vline, lab="",color=:black,line=(:dash, 1))
    xlabel!("Firing costs")
    ylabel!("Employment")
    xticks!([F1,F],[L"$\widehat{F}$",L"$F$"])
    xgrid!(:off)
    ygrid!(:off)
    savefig("ft_emp_color_"*f*".pdf")
    p3 = plot(grid_F, np,lab="Open-ended",color=:red,line=(:solid, 1))
    plot!(grid_F, u,lab="Non-employment",color=:black,line=(:solid, 1))
    plot!([F1,F], seriestype=:vline, lab="",color=:black,line=(:dash, 1))
    xlabel!("Firing costs")
    ylabel!("Employment")
    xticks!([F1,F],[L"$\widehat{F}$",L"$F$"])
    xgrid!(:off)
    ygrid!(:off)
    savefig("oe_emp_color_"*f*".pdf")

    # Productivity
    p4 = plot(grid_F, z_bar_p,lab=L"$\overline{z}^p$ - Dual",color=:red,line=(:solid, 1))
    plot!(grid_F, z_bar_f,lab=L"$\overline{z}^f$ - Dual",color=:blue,line=(:solid, 1))
    plot!(grid_F, z_bar,lab=L"$\overline{z}$ - Dual",color=:black,line=(:solid, 1))
    plot!(grid_F, z_bar_c,lab=L"$\overline{z}$ - Classic",color=:black,line=(:dash, 1))
    plot!([F1,F], seriestype=:vline, lab="",color=:black,line=(:dash, 1))
    xlabel!("Firing costs")
    ylabel!("Productivity")
    xticks!([F1,F],[L"$\widehat{F}$",L"$F$"])
    xgrid!(:off)
    ygrid!(:off)
    savefig("productivity_color_"*f*".pdf")

    # Welfare
    p5 = plot(grid_F, W,lab=L"Dual",color=:black,line=(:solid, 2))
    plot!(grid_F, W_c,lab=L"Classic",color=:black,line=(:dash, 2))
    plot!([F1,F], seriestype=:vline, lab="",color=:black,line=(:dash, 1))
    xlabel!("Firing costs")
    ylabel!("Welfare")
    xticks!([F1,F],[L"$\widehat{F}$",L"$F$"])
    xgrid!(:off)
    ygrid!(:off)
    savefig("welfare_color_"*f*".pdf")

    # Adjustment
    p6 = plot(grid_F, τ_p,lab="Open-ended",color=:red,line=(:solid, 1))
    plot!(grid_F, τ_f,lab="Fixed-term",color=:blue,line=(:solid, 1))
    # plot!(grid_F, τ_c,lab=L"$n$ - Classic",color=:black,line=(:dash, 2))
    plot!([F1,F], seriestype=:vline, lab="",color=:black,line=(:dash, 1))
    xlabel!("Firing costs")
    ylabel!(L"$\tau_{99\%}$")
    xticks!([F1,F],[L"$\widehat{F}$",L"$F$"])
    xgrid!(:off)
    ygrid!(:off)
    savefig("adj_colors_"*f*".pdf")

    cd("..")
end

"""
    dual_to_classic(F_ref, param, f; x0 = ... )

Plots the transition between a dual and a classic equilibrium when firing costs unexpectedly move from their value stored in `param` to `F_ref`. The resulting figures will have the suffix `f`.
"""
function dual_to_classic(F_ref, param, f; x0 = [0.9, 1.3])

    dual = DualSS()
    fill_DualSS!(dual, param, x0)
    param_ref = deepcopy(param)
    param_ref.F = F_ref
    classic = ClassicSS()
    fill_ClassicSS!(classic, param_ref, x0)

    @load_struct_cat old dual npᶜ np⁰ np nf zp zc zs u
    @load_struct classic nᶜ n⁰ ξ μ zp zc
    Macros.@load_struct param_ref δ α s λ
    @spec_funs

    @dual_to_classic
    cd("figures")

    t = range(0.,stop=24,length=1000)
    p1 = plot(t,vcat(np_old, nᵖ.(t[2:end])),lab=L"$n^p$",color=:black,line=(:solid, 1))
    plot!(t,vcat(nf_old, nᶠ.(t[2:end])),lab=L"$n^f$",color=:black,line=(:dot, 2))
    plot!(t,vcat(u_old, u.(t[2:end])),lab=L"$u$",color=:black,line=(:dash, 2), legend=:right)
    xlabel!("Time (months)")
    ylabel!("Employment")
    xgrid!(:off)
    ygrid!(:off)

    savefig("dual_to_classic_"*f*".pdf")

    p2 = plot(t,vcat(np_old, nᵖ.(t[2:end])),lab="Open-ended",color=:red,line=(:solid, 1))
    plot!(t,vcat(nf_old, nᶠ.(t[2:end])),lab="Fixed-term",color=:blue,line=(:dot, 1))
    plot!(t,vcat(u_old, u.(t[2:end])),lab="Non-employment",color=:black,line=(:dash, 1), legend=:right)
    xlabel!("Time (months)")
    ylabel!("Employment")
    xgrid!(:off)
    ygrid!(:off)

    savefig("dual_to_classic_color_"*f*".pdf")

    cd("..")

    return p1

end

"""
    classic_to_dual(F_ref, param, f; x0 = ... )

Plots the transition between a classic and a dual equilibrium when firing costs unexpectedly move from their value stored in `param` to `F_ref`. The resulting figures will have the suffix `f`.
"""
function classic_to_dual(param, f; x0 = [0.9, 1.3])

    classic = ClassicSS()
    fill_ClassicSS!(classic, param, x0)
    dual = DualSS()
    fill_DualSS!(dual, param, x0)

    @load_struct_cat old classic nᶜ n⁰ n zp zc
    npᶜ_old = nᶜ_old
    np⁰_old = n⁰_old
    np_old = n_old
    nf_old = 0.
    u_old = 1-np_old-nf_old
    @load_struct dual npᶜ np⁰ nf zc zs zf zp np ξ μᶠ μᵖ
    Macros.@load_struct param δ λ s α
    @spec_funs

    @classic_to_dual

    cd("figures")

    t = range(0,stop=60,length=1000)
    p1 = plot(t,vcat(np_old,nᵖ.(t[2:end])),lab=L"$n^p$",color=:black,line=(:solid, 1))
    plot!(t,vcat(nf_old, nᶠ.(t[2:end])),lab=L"$n^f$",color=:black,line=(:dot, 2))
    plot!(t,vcat(u_old,u.(t[2:end])),lab=L"$u$",color=:black,line=(:dash, 2), legend=:right)
    xlabel!("Time (months)")
    xgrid!(:off)
    ygrid!(:off)

    savefig("classic_to_dual_"*f*".pdf")

    cd("..")

    return p1
end

function plot_transitions(F_ref, param, f; x0 = [0.9, 1.3])
    p1 = dual_to_classic(F_ref, param, f; x0 = x0)
    p2 = classic_to_dual(param,f; x0 = x0)
    plot(p1,p2,layout=(1,2))
    cd("figures")
    savefig("dyn_reforms_"*f*".pdf")
    cd("..")
end

end


text:
struct BareCoefficient{T<:Integer}
    #=
    Symbolic representation of a coefficent in the bare Majorana Hamiltonian.
    =#

    site::T # the site of the corresponding bare operator
    type::String # the type of the corresponding bare operator

end

# Formatted printing
show(io::IO, bc::BareCoefficient) = print(io, bc.type, "_{", bc.site, "}")

# Needed for broadcasting mutliplication by bare coefficients
length(bc::BareCoefficient) = 1
iterate(bc::BareCoefficient) = (bc, 1)
iterate(bc::BareCoefficient, state::Int) = nothing

function isless(bc1::BareCoefficient, bc2::BareCoefficient)
    #=
    Are two bare coefficients ordered?
    =#

    # sort by site first, then by type
    bc1.site < bc2.site || (bc1.site == bc2.site && bc1.type < bc2.type)

end

function (==)(bc1::BareCoefficient, bc2::BareCoefficient)
    #=
    Do two bare coefficients represent the same value?
    =#

    # site and type must both match
    bc1.site == bc2.site && bc1.type == bc2.type

end

max_stars_repo_name: aormorningstar/SymbolicMajoranas

max_issues_repo_name: aormorningstar/SymbolicMajoranas

max_forks_repo_name: aormorningstar/SymbolicMajoranas


text:
include("load.jl")
StaticStoragesTests.runtests()


text:
module RoboPuzzle

function __init__()
    # silence the old "QApplication: invalid style override passed, ignoring
    #  it." warning...
    Base.eval(Main, :(ENV["QT_STYLE_OVERRIDE"] = ""))
end

using PyCall, PyPlot, Printf, LinearAlgebra

include("./Acorn/Acorn.jl")

include("./rbzl_constants.jl")
include("./rbzl_player.jl")
include("./rbzl_board.jl")
include("./rbzl_factory.jl")
include("./rbzl_helpers.jl")
include("./rbzl_acorn.jl")

function run(k::Integer=1)
    PyDict(pyimport("matplotlib")."rcParams")["toolbar"] = "None"
    load_level(k)
    idir = joinpath(@__DIR__, "..", "programs")
    if !isdir(idir)
        mkdir(idir)
    end
    ifile = joinpath(idir, @sprintf("level_%02d.jl", k))
    acorn_init()
    Acorn.acorn(ifile)
    close("all")
end

# ============================================================================ #
# GLOBALS
const BOARD = Vector{Board}(undef, 1)
const RAW_STR = Vector{String}(undef, 1)
# ============================================================================ #
end


text:
module DiffEqJumpExtensions

using DiffEqBase, Compat, Requires
using DiffEqJump

import DiffEqBase: DiscreteCallback, init, solve, solve!, plot_indices
import DiffEqJump: aggregate, get_jump_info_tuples, get_jump_info_fwrappers, executerx!,
                   evalrxrate, fill_cur_rates, register_next_jump_time!, build_jump_aggregation,
                   get_num_majumps, make_dependency_graph
import Base: size, getindex, setindex!, length, similar, indices, show, isempty, start, done, next
import Base.Order: Forward, Ordering, lt

include("aggregators/aggregators.jl")
include("aggregators/priorityqueue.jl")
#include("aggregators/simplepq.jl")
include("aggregators/nrm.jl")

export NRM, aggregate

end # module

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
# module Ecto

module Repo

# Repo.Schema.t
# Ecto.Schema.t
import ....Ecto

# Query.t
import ..Query

# Changeset.t
import ..Changeset

adapter = nothing

type t
    repo::Base.Random.UUID
    get::Function
    get_by::Function
    insert!::Function
    update!::Function
    insert_or_update!::Function
    delete!::Function
    update_all::Function
    delete_all::Function
    insert::Function
    update::Function
    insert_or_update::Function
    delete::Function
    function t()
        repo = Base.Random.uuid1()
        f(repo, g) = (args...; kw...) -> g(repo, args...; kw...)
        new(repo,
            f(repo, Repo.get),
            f(repo, Repo.get_by),
            f(repo, Repo.insert!),
            f(repo, Repo.update!),
            f(repo, Repo.insert_or_update!),
            f(repo, Repo.delete!),
            f(repo, Repo.update_all),
            f(repo, Repo.delete_all),
            f(repo, Repo.insert),
            f(repo, Repo.update),
            f(repo, Repo.insert_or_update),
            f(repo, Repo.delete)
        )
    end
end

include("repo/schema.jl")
include("repo/queryable.jl")

## Repo.Queryable
function get(repo::Base.Random.UUID, modul::Any, id::Any)#::Ecto.Schema.t
    opts = Dict()
    Repo.Queryable.get(repo, adapter, modul, id, opts)
end

function get_by(repo::Base.Random.UUID, modul::Any; kw...)#::Ecto.Schema.t
    opts = Dict(kw)
    get_by(repo, modul, opts)
end

function get_by(repo::Base.Random.UUID, modul::Any, opts::Dict)#::Ecto.Schema.t
    clause = nothing
    Repo.Queryable.get_by(repo, adapter, modul, clause, opts)
end

## Repo.Schema
function insert!(repo::Base.Random.UUID, struct_or_changeset::Union{Ecto.Schema.t,Changeset.t}; kw...)::Ecto.Schema.t
    Repo.Schema.insert!(repo, adapter, struct_or_changeset, Dict(kw))
end

function update!(repo::Base.Random.UUID, struct_or_changeset::Union{Ecto.Schema.t,Changeset.t}; kw...)::Ecto.Schema.t
    Repo.Schema.update!(repo, adapter, struct_or_changeset, Dict(kw))
end

function insert_or_update!(repo::Base.Random.UUID, struct_or_changeset::Union{Ecto.Schema.t,Changeset.t}; kw...)::Ecto.Schema.t
    Repo.Schema.insert_or_update!(repo, adapter, struct_or_changeset, Dict(kw))
end

function delete!(repo::Base.Random.UUID, struct_or_changeset::Union{Ecto.Schema.t,Changeset.t}; kw...)::Ecto.Schema.t
    Repo.Schema.delete!(repo, adapter, struct_or_changeset, Dict(kw))
end

## update_all
function update_all(repo::Base.Random.UUID, modul::Module; kw...)
# set=
end

function update_all(repo::Base.Random.UUID, modul::Module, updates::Array; kw...)
    opts = Dict(kw)
    if haskey(opts, :returning)
        isempty(opts[:returning]) && throw(ArgumentError("returning expects at least one field to be given"))
    end
end

function update_all{Q<:Query.Queryable}(repo::Base.Random.UUID, queryable::Q, updates::Array)
end

function update_all{Q<:Query.Queryable}(repo::Base.Random.UUID, queryable::Q; kw...)
    throw(Ecto.QueryError(""))
end

function update_all(repo::Base.Random.UUID, updates::Array)::Query.Applicable
    Query.Applicable()
end

function update_all(repo::Base.Random.UUID; kw...)::Query.Applicable
    Query.Applicable()
# set=
end

## delete_all
function delete_all(repo::Base.Random.UUID, modul::Module; kw...)
    opts = Dict(kw)
    if haskey(opts, :returning)
        isempty(opts[:returning]) && throw(ArgumentError("returning expects at least one field to be given"))
    end
end

function delete_all{Q<:Query.Queryable}(repo::Base.Random.UUID, queryable::Q)
    if !isempty(queryable, :select) || !isempty(queryable, :order_bys)
        throw(Ecto.QueryError(""))
    end
end

function insert(repo::Base.Random.UUID, changeset::Changeset.t)::Tuple{Symbol,Union{Changeset.t,Ecto.Schema.t}}
    opts = Dict()
    Repo.Schema.insert(repo, adapter, changeset, opts)
end

function update(repo::Base.Random.UUID, changeset::Changeset.t)::Tuple{Symbol,Union{Changeset.t,Ecto.Schema.t}}
    opts = Dict()
    Repo.Schema.update(repo, adapter, changeset, opts)
end

function insert_or_update(repo::Base.Random.UUID, changeset::Changeset.t)::Tuple{Symbol,Union{Changeset.t,Ecto.Schema.t}}
    opts = Dict()
    Repo.Schema.insert_or_update(repo, adapter, changeset, opts)
end

function delete(repo::Base.Random.UUID, changeset::Changeset.t)::Tuple{Symbol,Union{Changeset.t,Ecto.Schema.t}}
    opts = Dict()
    Repo.Schema.delete(repo, adapter, changeset, opts)
end

end # module Ecto.Repo


text:
# License for module ModiaMath.Logging: MIT
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module ModiaMath.Utilities

Utility functions used in ModiaMath.

# Main developer
Martin Otter, [DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module Utilities

export printobj

"""
    str = displayAsString(x)

Return the result of the Julia built-in command display(x) as string
instead of displaying it on the terminal.
"""
displayAsString(x) = sprint((io, x) -> display(TextDisplay(io), x), x)

"""
    str = displayAsString2(x)

Return the result of the Julia built-in command display(x) as string
instead of displaying it on the terminal, remove the preceding type
information and shift the whole output to the right
"""
function displayAsString2(obj)
    str  = displayAsString(obj)
    str2 = str[searchindex(str, "\n") + 1:end]
    str3 = "     " * replace(str2, "\n", "\n     ")
end

"""
    ModiaMath.printobj(name,obj)

Pretty print obj as "<name> = <display(obj)>".
"""
function printobj(name, obj)
    str = displayAsString2(obj)
    println("\n  ", name, " =\n", str)
end

end


text:
@testset "Bit buffer" begin
    bit_buffer = Tracking.BitBuffer()
    @test @inferred(get_bits(bit_buffer)) == 0
    @test @inferred(Tracking.length(bit_buffer)) == 0
    @test @allocated(Tracking.BitBuffer()) == 0

    prompt_accumulator = 0.0 + 0.0im
    secondary_code_or_bit_found = false
    prev_code_phase = 0.0
    code_phase = 1023
    prompt_correlator = 1.0 + 0.0im
    integration_time = 5ms
    gpsl1 = GPSL1()
    next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer(
        gpsl1,
        bit_buffer,
        prompt_accumulator,
        secondary_code_or_bit_found,
        prev_code_phase,
        code_phase,
        integration_time,
        prompt_correlator
    )

    @test next_prompt_accumulator == 0.0 + 0.0im
    @test @inferred(get_bits(next_bit_buffer)) == 0
    @test @inferred(Tracking.length(next_bit_buffer)) == 0

    prompt_accumulator = 20.0 + 0.0im
    secondary_code_or_bit_found = true
    prev_code_phase = 0.0
    code_phase = 1023
    prompt_correlator = 1.0 + 0.0im
    integration_time = 5ms

    next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer(
        gpsl1,
        bit_buffer,
        prompt_accumulator,
        secondary_code_or_bit_found,
        prev_code_phase,
        code_phase,
        integration_time,
        prompt_correlator
    )
    @test next_prompt_accumulator == 21.0 + 0.0im
    @test @inferred(get_bits(next_bit_buffer)) == 0
    @test @inferred(Tracking.length(next_bit_buffer)) == 0

    prompt_accumulator = 20.0 + 0.0im
    secondary_code_or_bit_found = true
    prev_code_phase = 9000
    code_phase = 1023
    prompt_correlator = 1.0 + 0.0im
    integration_time = 20ms

    next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer(
        gpsl1,
        bit_buffer,
        prompt_accumulator,
        secondary_code_or_bit_found,
        prev_code_phase,
        code_phase,
        integration_time,
        prompt_correlator
    )
    @test next_prompt_accumulator == 0.0 + 0.0im
    @test @inferred(get_bits(next_bit_buffer)) == 1
    @test @inferred(Tracking.length(next_bit_buffer)) == 1

    prompt_accumulator = 20.0 + 0.0im
    secondary_code_or_bit_found = true
    prev_code_phase = 0
    code_phase = 0
    prompt_correlator = 1.0 + 0.0im
    integration_time = 20ms

    next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer(
        gpsl1,
        bit_buffer,
        prompt_accumulator,
        secondary_code_or_bit_found,
        prev_code_phase,
        code_phase,
        integration_time,
        prompt_correlator
    )
    @test next_prompt_accumulator == 0.0 + 0.0im
    @test @inferred(get_bits(next_bit_buffer)) == 1
    @test @inferred(Tracking.length(next_bit_buffer)) == 1

end


text:
using StatsBase
import POMDPs: MDP, Policy

#
# Define a policy that will act as an crossentropy policy
#
type ValueTablePolicy <: Policy
    value_map::Array{Float64,2}
    action_space::Vector{Int64}
    epsilon::Float64
    learning_rate::Float64
    discount::Float64
end

#
# Default uniform intializer based on a games input
#
function ValueTablePolicy(pomdp::MDP, epsilon::Float64, learning_rate::Float64, discount::Float64)
    action_count = n_actions(pomdp)
    default_value_map = zeros(Float64, action_count, n_states(pomdp))
    return ValueTablePolicy(default_value_map, 1:action_count, epsilon, learning_rate, discount)
end

#
# Select random action from a policy based on a state and action weight
#
function action(policy::ValueTablePolicy, s::Int64)

    value, index = findmax(policy.value_map[:, s])
    is_uniform = iszero(policy.value_map[:, s])

    return if policy.epsilon < rand() && ~is_uniform policy.action_space[index] else rand(policy.action_space) end
end

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
module CommunalHelperCassetteSwapBlock

using ..Ahorn, Maple
using Ahorn.CommunalHelper

function swapFinalizer(entity)
    x, y = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    entity.data["nodes"] = [(x + width, y)]
end

@mapdef Entity "CommunalHelper/CassetteSwapBlock" CassetteSwapBlock(
    x::Integer,
    y::Integer,
    width::Integer=Maple.defaultBlockWidth,
    height::Integer=Maple.defaultBlockHeight,
    index::Integer=0,
    tempo::Number=1.0,
    noReturn::Bool=false,
    customColor="",
)

const ropeColors = Dict{Int,Ahorn.colorTupleType}(
    1 => (194, 116, 171, 255) ./ 255,
    2 => (227, 214, 148, 255) ./ 255,
    3 => (128, 224, 141, 255) ./ 255,
)

const defaultRopeColor = (110, 189, 245, 255) ./ 255

const placements = Ahorn.PlacementDict(
    "Cassette Swap Block ($index - $color) (Communal Helper)" => Ahorn.EntityPlacement(
        CassetteSwapBlock,
        "rectangle",
        Dict{String,Any}(
            "index" => index,
        ),
        swapFinalizer,
    ) for (color, index) in cassetteColorNames
)

Ahorn.editingOptions(entity::CassetteSwapBlock) = Dict{String,Any}(
    "index" => cassetteColorNames,
)

Ahorn.nodeLimits(entity::Maple.SwapBlock) = 1, 1

Ahorn.minimumSize(entity::CassetteSwapBlock) = 16, 16
Ahorn.resizable(entity::CassetteSwapBlock) = true, true

function Ahorn.selection(entity::CassetteSwapBlock)
    x, y = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return [
        Ahorn.Rectangle(x, y, width, height),
        Ahorn.Rectangle(stopX, stopY, width, height),
    ]
end

const block = "objects/cassetteblock/solid"
const crossSprite = "objects/CommunalHelper/cassetteMoveBlock/x"

function renderTrail(ctx, x::Number, y::Number, width::Number, height::Number, trail::String, index::Int, color=nothing)
    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    ropeColor = isnothing(color) ? get(ropeColors, index, defaultRopeColor) : color

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, trail, x + (i - 1) * 8, y + 2, 6, 0, 8, 6, tint=ropeColor)
        Ahorn.drawImage(ctx, trail, x + (i - 1) * 8, y + height - 8, 6, 14, 8, 6, tint=ropeColor)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, trail, x + 2, y + (i - 1) * 8, 0, 6, 6, 8, tint=ropeColor)
        Ahorn.drawImage(ctx, trail, x + width - 8, y + (i - 1) * 8, 14, 6, 6, 8, tint=ropeColor)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, trail, x + (i - 1) * 8, y + (j - 1) * 8, 6, 6, 8, 8, tint=ropeColor)
    end

    Ahorn.drawImage(ctx, trail, x + width - 8, y + 2, 14, 0, 6, 6, tint=ropeColor)
    Ahorn.drawImage(ctx, trail, x + width - 8, y + height - 8, 14, 14, 6, 6, tint=ropeColor)
    Ahorn.drawImage(ctx, trail, x + 2, y + 2, 0, 0, 6, 6, tint=ropeColor)
    Ahorn.drawImage(ctx, trail, x + 2, y + height - 8, 0, 14, 6, 6, tint=ropeColor)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::CassetteSwapBlock, room::Maple.Room)
    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    index = Int(get(entity.data, "index", 0))

    hexColor = String(get(entity.data, "customColor", ""))
    if hexColor != "" && length(hexColor) == 6
        renderCassetteBlock(ctx, stopX, stopY, width, height, index, hexToRGBA(hexColor))
    else
        renderCassetteBlock(ctx, stopX, stopY, width, height, index)
    end

    Ahorn.drawArrow(ctx, startX + width / 2, startY + height / 2, stopX + width / 2, stopY + height / 2, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.renderAbs(
    ctx::Ahorn.Cairo.CairoContext,
    entity::CassetteSwapBlock,
    room::Maple.Room,
)
    startX, startY = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    index = Int(get(entity.data, "index", 0))
    color = getCassetteColor(index)

    hexColor = String(get(entity.data, "customColor", ""))
    if hexColor != "" && length(hexColor) == 6
        color = hexToRGBA(hexColor)
    end

    renderTrail(ctx, min(startX, stopX), min(startY, stopY), abs(startX - stopX) + width, abs(startY - stopY) + height, "objects/swapblock/target", index, color)

    renderCassetteBlock(ctx, startX, startY, width, height, index, color)

    if Bool(get(entity.data, "noReturn", false))
        noReturnSprite = Ahorn.getSprite(crossSprite, "Gameplay")
        Ahorn.drawImage(ctx, noReturnSprite, startX + div(width - noReturnSprite.width, 2), startY + div(height - noReturnSprite.height, 2), tint=color)
    end
end

end

max_stars_repo_name: vitellaryjr/CommunalHelper

max_issues_repo_name: vitellaryjr/CommunalHelper

max_forks_repo_name: vitellaryjr/CommunalHelper


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

for i in [(:_count_a, :a_bitcount), (:_count_c, :c_bitcount), (:_count_g, :g_bitcount), (:_count_t, :t_bitcount)]
    @eval begin
        @inline function $(i[1])(alph::A, head::UInt64, tail...) where {A<:NucleicAcidAlphabet}
            return $(i[2])(head, alph) + $(i[1])(alph, tail...)
        end
        @inline $(i[1])(alph::A) where {A<:NucleicAcidAlphabet} = 0
    end
end

@inline function _count_gc(alph::A, head::UInt64, tail...) where {A<:NucleicAcidAlphabet}
    return gc_bitcount(head, alph) + _count_gc(alph, tail...)
end
@inline _count_gc(::A) where {A<:NucleicAcidAlphabet} = 0

count_a(x::Kmer{A,K,N})  where {A<:NucleicAcidAlphabet,K,N} = _count_a(A(), x.data...) - n_unused(x)
count_c(x::Kmer{A,K,N})  where {A<:NucleicAcidAlphabet,K,N} = _count_c(A(), x.data...)
count_g(x::Kmer{A,K,N})  where {A<:NucleicAcidAlphabet,K,N} = _count_g(A(), x.data...)
count_t(x::Kmer{A,K,N})  where {A<:NucleicAcidAlphabet,K,N} = _count_t(A(), x.data...)

count_gc(x::Kmer{A,K,N}) where {A<:NucleicAcidAlphabet,K,N} = _count_gc(A(), x.data...)
Base.count(::typeof(isGC), x::Kmer{A,K,N}) where {A<:NucleicAcidAlphabet,K,N} = count_gc(x)

# TODO: Expand to Amino Acid Kmers as well...
@inline function Base.count(::typeof(!=), a::Kmer{A,K,N}, b::Kmer{A,K,N}) where {A<:NucleicAcidAlphabet,K,N}
    ad = a.data
    bd = b.data
    sum = 0
    @inbounds for i in 1:N
        sum += BioSequences.mismatch_bitcount(ad[i], bd[i], A())
    end
    return sum
end

# TODO: Expand to Amino Acid Kmers as well...
@inline function Base.count(::typeof(==), a::Kmer{A,K,N}, b::Kmer{A,K,N}) where {A<:NucleicAcidAlphabet,K,N}
    ad = a.data
    bd = b.data
    sum = 0
    @inbounds for i in 1:N
        sum += BioSequences.match_bitcount(ad[i], bd[i], A())
    end
    return sum - n_unused(a)
end


text:
function simplify(c::Chain; rules = naive, mod = @__MODULE__)
	eg = EGraph(c)

	settermtype!(eg, :Chain, 0, Chain)
	settermtype!(eg, :Chain, 1, Chain)
	settermtype!(eg, :Chain, 2, Chain)
	settermtype!(eg, :Gate, 2, Gate)
	settermtype!(eg, :Ctrl, 2, Ctrl)

	saturate!(eg, rules; mod = @__MODULE__)
	c_op = extract!(eg, astsize)
	return postprocess(c_op)
end


text:
@with_kw mutable struct Iterate{T}
    x::Vector{T}
    y::Vector{T}
    z::Vector{T}
    s::Vector{T}
    τ::T
    κ::T
end

# τ0 = κ0 = 1 specified in CVXOPT User Guide, §7.1, pg 17:
Iterate{T}(x, y, z, s) where T = Iterate{T}(x, y, z, s, 1., 1.,)

function Iterate(prob::Problem{T}) where T
    x = zeros(T, get_num_vars(prob))
    y = zeros(T, get_num_constrs(prob))
    s = get_id(prob.C)
    z = get_id(prob.C)
    return Iterate{T}(x, y, z, s)
end

copy(iter::Iterate) = Iterate(copy(iter.x), copy(iter.y), copy(iter.z), copy(iter.s), copy(iter.τ), copy(iter.κ))

function copy_to!(iter::Iterate, best_iter::Iterate)
    best_iter.x = iter.x
    best_iter.y = iter.y
    best_iter.z = iter.z
    best_iter.s = iter.s
    best_iter.τ = iter.τ
    best_iter.κ = iter.κ
end

function is_valid(iter::Iterate{T}, prob::Problem{T}) where T
    return (all(isfinite.(iter.x)) && all(isfinite.(iter.y)) && all(isfinite.(iter.z)) &&
            all(isfinite.(iter.s)) && isfinite(iter.τ) && isfinite(iter.κ) &&
            iter.s ∈ prob.C && iter.z ∈ prob.C)
end

@with_kw struct Residual{T}
    rx::Vector{T}
    ry::Vector{T}
    rz::Vector{T}
    rτ::T

    @assert !any(isnan.(rx))
    @assert !any(isnan.(ry))
    @assert !any(isnan.(rz))
    @assert !any(isnan.(rτ))
end

stack(resid::Residual) = [resid.rx, resid.ry, resid.rz, resid.rτ]

@with_kw struct StepDirection{T}
    Δx::Vector{T}
    Δy::Vector{T}
    Δz::Vector{T}
    Δτ::T
    Δs::Vector{T}
    Δκ::T

    #@assert !any(isnan.(Δx))
    #@assert !any(isnan.(Δy))
    #@assert !any(isnan.(Δz))
    #@assert !any(isnan.(Δτ))
    #@assert !any(isnan.(Δs))
    #@assert !any(isnan.(Δκ))
    #@assert length(Δs) == length(Δz)
end

function StepDirection(Δx::Vector{T}, Δy::Vector{T}, Δz::Vector{T}, Δτ::Vector{T},
                       Δs::Vector{T}, Δκ::Vector{T}) where {T <: AbstractFloat}
    return StepDirection(Δx, Δy, Δz, Δτ[1], Δs, Δκ[1])
end

is_valid(dir::StepDirection) = (all(isfinite.(dir.Δx)) && all(isfinite.(dir.Δy)) && all(isfinite.(dir.Δz)) &&
                                all(isfinite.(dir.Δs)) && isfinite(dir.Δτ) && isfinite(dir.Δκ))


text:
# simple Julia APIs over Kubernetes Swagger interface

function sel(label::String, op::Symbol)
    @assert op === :exists
    label
end
sel(label::String, op::Symbol, items::String...) = label * " " * string(op) * " (" * join(items, ",") * ")"
sel(cnd::String...) = join(cnd, ", ")

_delopts(; kwargs...) = Typedefs.MetaV1.DeleteOptions(; preconditions=Typedefs.MetaV1.Preconditions(; kwargs...), kwargs...)

_kubectx(ctx::KuberContext) = ctx
_kubectx(ctx::KuberWatchContext) = ctx.ctx

function _get_apictx(ctx::Union{KuberContext,KuberWatchContext}, O::Symbol, apiversion::Union{String,Nothing}; max_tries::Int=retries(ctx, false))
    kubectx = _kubectx(ctx)
    kubectx.initialized || set_api_versions!(kubectx; max_tries=max_tries)

    if apiversion !== nothing
        k = Kuber.api_group_type(ctx, apiversion)
        apictx = k(kubectx.client)
    else
        kapi = kubectx.modelapi[O]
        apictx = kapi.api(kubectx.client)
    end
    apictx
end

_api_function(ctx::Union{KuberContext,KuberWatchContext}, name::Symbol) = isdefined(apimodule(ctx), name) ? apimodule(ctx).eval(name) : nothing
_api_function(ctx::Union{KuberContext,KuberWatchContext}, name) = _api_function(ctx, Symbol(name))

function watch(fn::Function, ctx::KuberContext; buffersize::Int=1024, stream::KuberEventStream=KuberEventStream(buffersize))
    watchctx = KuberWatchContext(ctx, stream)
    fn(watchctx, stream)
end

function watch(streamprocessor::Function, ctx::KuberContext, watched::Function, args...; kwargs...)
    watch(ctx) do watchctx, stream
        @sync begin
            @async try
                watched(watchctx, args...; kwargs...)
            finally
                close(stream)
            end
            @async streamprocessor(stream)
        end
    end
end

function list(ctx::Union{KuberContext,KuberWatchContext}, O::Symbol, name::String;
        apiversion::Union{String,Nothing}=nothing,
        namespace::Union{String,Nothing}=_kubectx(ctx).namespace,
        max_tries::Int=retries(ctx, false),
        watch=isa(ctx, KuberWatchContext),
        resourceVersion=nothing,
        kwargs...)
    apictx = _get_apictx(ctx, O, apiversion; max_tries=max_tries)
    namespaced = (namespace !== nothing) && !isempty(namespace)
    allnamespaces = namespaced && (namespace == "*")

    eventstream = isa(ctx, KuberWatchContext) ? ctx.stream : nothing
    result = nothing
    args = Any[name]
    if allnamespaces
        apicall = apimodule(ctx).eval(Symbol("list$(O)ForAllNamespaces"))
    elseif namespaced
        apicall = apimodule(ctx).eval(Symbol("listNamespaced$O"))
        push!(args, namespace)
    else
        apicall = apimodule(ctx).eval(Symbol("list$O"))
    end

    if !watch || resourceVersion === nothing
        result = k8s_retry(; max_tries=max_tries) do
            apicall(apictx, args...; kwargs...)
        end
    end

    # if not watching, return the first result
    watch || (return result)
    if result !== nothing
        resourceVersion = result.metadata.resourceVersion
        # push the first Event consisting of existing data
        put!(eventstream, result)
    end

    # start watch and return the HTTP response object on completion
    return k8s_retry(; max_tries=max_tries) do
        apicall(apictx, eventstream, args...; watch=watch, resourceVersion=resourceVersion, kwargs...)
    end
end

function list(ctx::Union{KuberContext,KuberWatchContext}, O::Symbol;
        apiversion::Union{String,Nothing}=nothing,
        namespace::Union{String,Nothing}=_kubectx(ctx).namespace,
        max_tries::Int=retries(ctx, false),
        watch=isa(ctx, KuberWatchContext),
        resourceVersion=nothing,
        kwargs...)
    apictx = _get_apictx(ctx, O, apiversion; max_tries=max_tries)
    namespaced = (namespace !== nothing) && !isempty(namespace)
    allnamespaces = namespaced && (namespace == "*")

    eventstream = isa(ctx, KuberWatchContext) ? ctx.stream : nothing

    result = nothing
    args = Any[]
    if allnamespaces
        apicall = apimodule(ctx).eval(Symbol("list$(O)ForAllNamespaces"))
    elseif namespaced
        apicall = apimodule(ctx).eval(Symbol("listNamespaced$O"))
        push!(args, namespace)
    else
        apicall = apimodule(ctx).eval(Symbol("list$O"))
    end

    if !watch || resourceVersion === nothing
        result = k8s_retry(; max_tries=max_tries) do
            apicall(apictx, args...; kwargs...)
        end
    end

    # if not watching, retuen the first result
    watch || (return result)
    if result !== nothing
        resourceVersion = result.metadata.resourceVersion
        # push the first Event consisting of existing data
        put!(eventstream, result)
    end

    # start watch and return the HTTP response object on completion
    return k8s_retry(; max_tries=max_tries) do
        apicall(apictx, eventstream, args...; watch=watch, resourceVersion=resourceVersion, kwargs...)
    end
end

function get(ctx::Union{KuberContext,KuberWatchContext}, O::Symbol, name::String;
        apiversion::Union{String,Nothing}=nothing,
        max_tries::Integer=retries(ctx),
        watch=isa(ctx, KuberWatchContext),
        resourceVersion=nothing,
        kwargs...)
    apictx = _get_apictx(ctx, O, apiversion; max_tries=max_tries)

    eventstream = isa(ctx, KuberWatchContext) ? ctx.stream : nothing
    result = nothing
    args = Any[name]
    if (apicall = _api_function(ctx, "read$O")) !== nothing
        # nothing
    elseif (apicall = _api_function(ctx, "readNamespaced$O")) !== nothing
        push!(args, _kubectx(ctx).namespace)
    else
        throw(ArgumentError("No API functions could be located using :$O"))
    end

    if !watch || resourceVersion === nothing
        result = k8s_retry(; max_tries=max_tries) do
            apicall(apictx, args...; kwargs...)
        end
    end

    # if not watching, retuen the first result
    watch || (return result)
    if result !== nothing
        resourceVersion = result.metadata.resourceVersion
        # push the first Event consisting of existing data
        put!(eventstream, result)
    end

    # start watch and return the HTTP response object on completion
    return k8s_retry(; max_tries=max_tries) do
        apicall(apictx, eventstream, args...; watch=watch, resourceVersion=resourceVersion, kwargs...)
    end
end

function get(ctx::Union{KuberContext,KuberWatchContext}, O::Symbol;
        apiversion::Union{String,Nothing}=nothing,
        label_selector=nothing,
        namespace::Union{String,Nothing}=_kubectx(ctx).namespace,
        max_tries::Integer=retries(ctx, false),
        watch=isa(ctx, KuberWatchContext),
        resourceVersion=nothing,
        kwargs...)
    apictx = _get_apictx(ctx, O, apiversion; max_tries=max_tries)

    eventstream = isa(ctx, KuberWatchContext) ? ctx.stream : nothing
    result = nothing
    args = Any[]
    apiname = "list$O"
    namespace === nothing && (apiname *= "ForAllNamespaces")
    if (apicall = _api_function(ctx, apiname)) !== nothing
        # nothing
    elseif (apicall = _api_function(ctx, "listNamespaced$O")) !== nothing
        push!(args, namespace)
    else
        throw(ArgumentError("No API functions could be located using :$O"))
    end

    if !watch || resourceVersion === nothing
        result = k8s_retry(; max_tries=max_tries) do
            apicall(apictx, args...; labelSelector=label_selector, kwargs...)
        end
    end

    # if not watching, retuen the first result
    watch || (return result)
    if result !== nothing
        resourceVersion = result.metadata.resourceVersion
        # push the first Event consisting of existing data
        put!(eventstream, result)
    end

    # start watch and return the HTTP response object on completion
    return k8s_retry(; max_tries=max_tries) do
        apicall(apictx, eventstream, args...; watch=watch, resourceVersion=resourceVersion, labelSelector=label_selector, kwargs...)
    end
end

function watch(ctx::KuberContext, O::Symbol, outstream::Channel, name::String;
        apiversion::Union{String,Nothing}=nothing,
        namespace::Union{String,Nothing}=ctx.namespace,
        max_tries::Int=retries(ctx, false),
        kwargs...)
    apictx = _get_apictx(ctx, O, apiversion; max_tries=max_tries)
    namespaced = (namespace !== nothing) && !isempty(namespace)
    allnamespaces = namespaced && (namespace == "*")

    if allnamespaces
        apicall = apimodule(ctx).eval(Symbol("watch$(O)ForAllNamespaces"))
        return k8s_retry(; max_tries=max_tries) do
            apicall(apictx, outstream, name; kwargs...)
        end
    elseif namespaced
        apicall = apimodule(ctx).eval(Symbol("watchNamespaced$O"))
        return k8s_retry(; max_tries=max_tries) do
            apicall(apictx, outstream, name, namespace; kwargs...)
        end
    else
        apicall = apimodule(ctx).eval(Symbol("watch$O"))
        return k8s_retry(; max_tries=max_tries) do
            apicall(apictx, outstream, name; kwargs...)
        end
    end
end

function watch(ctx::KuberContext, O::Symbol, outstream::Channel;
        apiversion::Union{String,Nothing}=nothing,
        namespace::Union{String,Nothing}=ctx.namespace,
        max_tries::Int=retries(ctx, false),
        kwargs...)
    apictx = _get_apictx(ctx, O, apiversion; max_tries=max_tries)
    namespaced = (namespace !== nothing) && !isempty(namespace)
    allnamespaces = namespaced && (namespace == "*")

    if allnamespaces
        apicall = apimodule(ctx).eval(Symbol("watch$(O)ForAllNamespaces"))
        return k8s_retry(; max_tries=max_tries) do
            apicall(apictx, outstream; kwargs...)
        end
    elseif namespaced
        apicall = apimodule(ctx).eval(Symbol("watchNamespaced$O"))
        return k8s_retry(; max_tries=max_tries) do
            apicall(apictx, outstream, namespace; kwargs...)
        end
    else
        apicall = apimodule(ctx).eval(Symbol("watch$O"))
        return k8s_retry(; max_tries=max_tries) do
            apicall(apictx, outstream; kwargs...)
        end
    end
end

function put!(ctx::KuberContext, v::T; max_tries::Int=retries(ctx, true)) where {T<:SwaggerModel}
    vjson = convert(Dict{String,Any}, v)
    put!(ctx, Symbol(vjson["kind"]), vjson; max_tries=max_tries)
end

function put!(ctx::KuberContext, O::Symbol, d::Dict{String,Any}; max_tries::Int=retries(ctx, true))
    apictx = _get_apictx(ctx, O, get(d, "apiVersion", nothing))
    if (apicall = _api_function(ctx, "create$O")) !== nothing
        return k8s_retry(; max_tries=max_tries) do
            apicall(apictx, d)
        end
    elseif (apicall = _api_function(ctx, "createNamespaced$O")) !== nothing
        return k8s_retry(; max_tries=max_tries) do
            apicall(apictx, ctx.namespace, d)
        end
    else
        throw(ArgumentError("No API functions could be located using :$O"))
    end
end

function delete!(ctx::KuberContext, v::T; max_tries::Int=retries(ctx, true), kwargs...) where {T<:SwaggerModel}
    vjson = convert(Dict{String,Any}, v)
    kind = vjson["kind"]
    name = vjson["metadata"]["name"]
    delete!(ctx, Symbol(kind), name; apiversion=get(vjson, "apiVersion", nothing), max_tries=max_tries, kwargs...)
end

function delete!(ctx::KuberContext, O::Symbol, name::String; apiversion::Union{String,Nothing}=nothing, max_tries::Int=retries(ctx, true), kwargs...)
    apictx = _get_apictx(ctx, O, apiversion)

    params = [apictx, name]

    if (apicall = _api_function(ctx, "delete$O")) !== nothing
        return k8s_retry(; max_tries=max_tries) do
            apicall(params...; kwargs...)
        end
    elseif (apicall = _api_function(ctx, "deleteNamespaced$O")) !== nothing
        push!(params, ctx.namespace)
        return k8s_retry(; max_tries=max_tries) do
            apicall(params...; kwargs...)
        end
    else
        throw(ArgumentError("No API functions could be located using :$O"))
    end
end

function update!(ctx::KuberContext, v::T, patch, patch_type; max_tries::Int=retries(ctx, true)) where {T<:SwaggerModel}
    vjson = convert(Dict{String,Any}, v)
    kind = vjson["kind"]
    name = vjson["metadata"]["name"]
    update!(ctx, Symbol(kind), name, patch, patch_type; apiversion=get(vjson, "apiVersion", nothing), max_tries=max_tries)
end

function update!(ctx::KuberContext, O::Symbol, name::String, patch, patch_type; apiversion::Union{String,Nothing}=nothing, max_tries::Int=retries(ctx, true))
    apictx = _get_apictx(ctx, O, apiversion)

    if (apicall = _api_function(ctx, "patch$O")) !== nothing
        return k8s_retry(; max_tries=max_tries) do
            apicall(apictx, name, patch; _mediaType=patch_type)
        end
    elseif (apicall = _api_function(ctx, "patchNamespaced$O")) !== nothing
        return k8s_retry(; max_tries=max_tries) do
            apicall(apictx, name, ctx.namespace, patch; _mediaType=patch_type)
        end
    else
        throw(ArgumentError("No API functions could be located using :$O"))
    end
end

"""
Fetch logs of the specified Pod

Parameters:
- name::String (required)

Keyword Args:
- container::String : The container for which to stream logs. Defaults to only container if there is one container in the pod.
- follow::Bool : Follow the log stream of the pod. Defaults to false.
- limitBytes::Int32 : If set, the number of bytes to read from the server before terminating the log output. This may not display a complete final line of logging, and may return slightly more or slightly less than the specified limit.
- pretty::String : If 'true', then the output is pretty printed.
- previous::Bool : Return previous terminated container logs. Defaults to false.
- sinceSeconds::Int32 : A relative time in seconds before the current time from which to show logs. If this value precedes the time a pod was started, only logs since the pod start will be returned. If this value is in the future, no logs will be returned. Only one of sinceSeconds or sinceTime may be specified.
- sinceTime::String : An RFC3339 timestamp from which to show logs. If this value precedes the time a pod was started, only logs since the pod start will be returned. If this value is in the future, no logs will be returned. Only one of sinceSeconds or sinceTime may be specified.
- tailLines::Int32 : If set, the number of lines from the end of the logs to show. If not specified, logs are shown from the creation of the container or sinceSeconds or sinceTime
- timestamps::Bool : If true, add an RFC3339 or RFC3339Nano timestamp at the beginning of every line of log output. Defaults to false.

Returns: String of all log entries, one per line
"""
get_logs(ctx::KuberContext, pod_name::String; kwargs...) = get(ctx, :PodLog, pod_name; kwargs...)

list_namespaced_custom_metrics(ctx::KuberContext, metricname::String; kwargs...) = list(ctx, :MetricValue, "metrics/"*metricname; kwargs...)
list_namespaced_custom_metrics(ctx::KuberContext, objecttype::String, metricname::String; kwargs...) = list(ctx, :MetricValue, objecttype * "/*/" * metricname; kwargs...)
list_namespaced_custom_metrics(ctx::KuberContext, objecttype::String, objectname::String, metricname::String; kwargs...) = list(ctx, :MetricValue, objecttype * "/" * objectname * "/" * metricname; kwargs...)

list_custom_metrics(ctx::KuberContext, objecttype::String, metricname::String; kwargs...) = list(ctx, :MetricValue, objecttype * "/*/" * metricname; namespace=nothing, kwargs...)
list_custom_metrics(ctx::KuberContext, objecttype::String, objectname::String, metricname::String; kwargs...) = list(ctx, :MetricValue, objecttype * "/" * objectname * "/" * metricname; namespace=nothing, kwargs...)


text:
@with_kw struct sparameters{R}
  #thermal constants
  b1::R = 0.243;
  b2::R = 0.393;
  b3::R = 1.534;
  #retention - VG model
  alpha::R = 5;
  n::R = 2;
  m::R = 1-1/n;
  theta_s::R = 0.4;
  theta_r::R = 0.01;
  Ks::R = 2/(24*60*60);
  #matric head top constrain
  h_max::R = -1e20;
  #additional parameters
  G_wT::R = 7; # gain factor [-] taken from Jiangbo et al. (2017)
  fc::R = 0.02; #taken from Jiangbo et al. (2017)
end

global soil_parameters=sparameters{Float64}()


text:
using Metatheory
using Metatheory.EGraphs
using Test

@metatheory_init ()

struct MyExpr
    head::Any
    # NOTE! this will not work, when replacing
    # with z in the theory defined below, the arg type
    # will be EGraphs.EClass! Additional manipulation
    # is needed for custom term types with stricter arg types
    # args::Vector{Union{Int, MyExpr}}
    args::Vector{Any}
    # additional metadata
    foo::String
    bar::Vector{Complex}
    baz::Set{Int}
end

import Base.(==)
(==)(a::MyExpr, b::MyExpr) = a.head == b.head && a.args == b.args &&
    a.foo == b.foo && a.bar == b.bar && a.baz == b.baz

MyExpr(head, args) = MyExpr(head, args, "", Complex[], Set{Int}())
MyExpr(head) = MyExpr(head, [])

# Methods needed by `src/TermInterface.jl`
TermInterface.gethead(e::MyExpr) = e.head
TermInterface.getargs(e::MyExpr) = e.args
TermInterface.istree(e::MyExpr) = true
# NamedTuple
TermInterface.getmetadata(e::MyExpr) = (foo=e.foo, bar=e.bar, baz=e.baz)
TermInterface.preprocess(e::MyExpr) = MyExpr(e.head, e.args, uppercase(e.foo), e.bar, e.baz)

# f(g(2), h(4)) with some metadata in h
hcall = MyExpr(:call, [:h, 4], "hello", [2+3im, 4+2im], Set{Int}([4,5,6]))
ex = MyExpr(:call, [:f, MyExpr(:call, [:g, 2]), hcall])

function EGraphs.instantiateterm(g::EGraph, pat::PatTerm,  T::Type{MyExpr}, sub::Sub, rule::Rule)
    head = pat.head
    args_start = 1
    # if pat.head == :call
    #     pl = pat.args[1]
    #     head = pl.val
    #     args_start = 2
    # end

    args = map(x -> EGraphs.instantiate(g, x, sub, rule), (@view pat.args[args_start:end]))

    res = MyExpr(head, args, "", Complex[], Set{Int64}())
    println(res)
    res
end

# Define an extraction method dispatching on MyExpr
function EGraphs.extractnode(n::ENode{MyExpr}, extractor::Function)
    (foo, bar, baz) = n.metadata
    # extracted arguments
    ret_args = []

    for a ∈ n.args
        push!(ret_args, extractor(a))
    end

    return MyExpr(n.head, ret_args, foo, bar, baz)
end

# let's create an egraph
g = EGraph(ex)

# let's create an example theory
t = @theory begin
    # this way, z will be a regular expr
    # f(g(2), a) => z(a)
    # we can use dynamic rules to construct values of type MyExpr
    # f(g(2), a) |> MyExpr(:z, [a])

    # terms in the RHS inherit the type of terms in the lhs
    f(g(2), a) => f(a)
end

saturate!(g, t; mod=@__MODULE__)

display(g.classes)

expected = MyExpr(:call, [:f, MyExpr(:call, [:h, 4], "HELLO", Complex[2 + 3im, 4 + 2im], Set([5, 4, 6]))], "", Complex[], Set{Int64}())

extracted = extract!(g, astsize)

@test expected == extracted


text:
using Revise

Revise.revise()

using Term
import Term: install_stacktrace

install_stacktrace()

println(@green "Metaprogramming is great!")
println(@underline "Made with [green]Term.jl[/green]")

println(
    @style "This is my text, effortlessly [black bold on_red]styled![/black bold on_red]" bold blue underline
)


text:
@testset "Hierarchies" begin

  # need expanding
  hierarchy_result = Hierarchy("Hierarchy 1")
  @test hierarchy_result.Name == "Hierarchy 1"
  @test name(hierarchy_result) == "Hierarchy 1"

  hierarchy_2 = JSON.parse("""
  {
    "Name": "Hierarchy 2"
  }
  """)

  @test name(Hierarchy(hierarchy_2)) == "Hierarchy 2"
  @test Hierarchy(hierarchy_2) != hierarchy_result

end


text:
# note that this function assumes that the number of nonzeros, number of cols, and number of rows is already the same
function Base.copyto!(dest::SuiteSparse.CHOLMOD.Sparse, src::SparseMatrixCSC)
	crep = unsafe_load(dest.ptr)
	ht,wd = size(src)
	nnz = length(src.nzval)
	@assert crep.nzmax == nnz "$(crep.nzmax) neq $(nnz)"
	@assert crep.nrow == src.m
	@assert crep.ncol == src.n
	for i=1:src.n+1
		unsafe_store!(crep.p, (@inbounds src.colptr[i])-1, i)
	end
	for i=1:nnz
		unsafe_store!(crep.i, (@inbounds src.rowval[i])-1, i)
		unsafe_store!(crep.x, (@inbounds src.nzval[i]), i)
	end
end

function Base.copyto!(dest::SuiteSparse.CHOLMOD.Dense{T}, src::AbstractArray{T, 1}) where {T}
	uref = unsafe_load(dest.ptr)
	@assert uref.nrow == length(src) "$(uref.nrow) neq $(length(src))"
	for i=1:length(src)
		unsafe_store!(uref.x, src[i], i)
	end
end

function LinearAlgebra.lmul!(l::Diagonal, r::SuiteSparse.CHOLMOD.Sparse)
	crep = unsafe_load(r.ptr)
	diag = l.diag
	for i=1:crep.ncol
		colstart = unsafe_load(crep.p, i)+1
		colend = unsafe_load(crep.p, i+1)
		for j=colstart:colend
			rowidx = unsafe_load(crep.i, j)
			cellv = unsafe_load(crep.x, j)
			unsafe_store!(crep.x, cellv*diag[i], j)
		end
	end
end

struct LdivWorkspace{T<:SuiteSparse.CHOLMOD.VTypes}
	Xref::Ref{Ptr{Nothing}}
	Y::Ref{Ptr{Nothing}}
	E::Ref{Ptr{Nothing}}
	LdivWorkspace(::Type{T}) where T = new{T}(Ref(C_NULL), Ref(C_NULL), Ref(C_NULL))
end

function div!(F::SuiteSparse.CHOLMOD.Factor{Tv}, x::SuiteSparse.CHOLMOD.Dense{Tv}, B::SuiteSparse.CHOLMOD.Dense{Tv}, W::LdivWorkspace{Tv}) where Tv<:SuiteSparse.CHOLMOD.VTypes
    if size(F,1) != size(B,1)
        throw(DimensionMismatch("LHS and RHS should have the same number of rows. " *
            "LHS has $(size(F,1)) rows, but RHS has $(size(B,1)) rows."))
    end
    if !issuccess(F)
        s = unsafe_load(pointer(F))
        if s.is_ll == 1
            throw(LinearAlgebra.PosDefException(s.minor))
        else
            throw(LinearAlgebra.ZeroPivotException(s.minor))
        end
    end
    W.Xref[] = Base.unsafe_convert(Ptr{Nothing}, Base.unsafe_convert(Ptr{SuiteSparse.CHOLMOD.C_Dense{Tv}}, x))
    res = ccall((SuiteSparse.CHOLMOD.@cholmod_name("solve2"),:libcholmod), Cint,
            (	# input
            	Cint, Ptr{SuiteSparse.CHOLMOD.C_Factor{Tv}}, Ptr{SuiteSparse.CHOLMOD.C_Dense{Tv}}, Ptr{Nothing},
            	# output
            	Ref{Ptr{Nothing}}, Ptr{Nothing},
            	# workspace
            	Ref{Ptr{Nothing}}, Ref{Ptr{Nothing}},
            	# common
            	Ptr{UInt8}),
                SuiteSparse.CHOLMOD.CHOLMOD_A, F, B, C_NULL,
                W.Xref, C_NULL,
                W.Y, W.E,
                SuiteSparse.CHOLMOD.common_struct[Threads.threadid()])
    return res
end

function change_xtype!(F::SuiteSparse.CHOLMOD.Sparse{T}, xtype::Int32) where T
	return ccall((SuiteSparse.CHOLMOD.@cholmod_name("sparse_xtype"), :libcholmod), Cint,
		(Cint, Ptr{SuiteSparse.CHOLMOD.C_Sparse{T}}, Ptr{UInt8}),
		xtype, F, SuiteSparse.CHOLMOD.common_struct[Threads.threadid()])
end


text:
using JSON
include("utils.jl")

function notify_overheads()
    results = [Ref{Any}() for _ in 1:Threads.nthreads()]
    counter1 = Ref(0)
    cond = Threads.Condition()
    tasks = Task[]
    for (tid, ref) in pairs(results)
        task = @task begin
            @assert tid == Threads.threadid()
            # print('.')
            lock(cond) do
                if (counter1[] += 1) == Threads.nthreads()
                    # println("notifying")
                    t0 = time_ns()
                    notify(cond)
                    t1 = time_ns()
                    ref[] = (t0, t1)
                else
                    while counter1[] < Threads.nthreads()
                        wait(cond)
                    end
                    ref[] = time_ns()
                end
            end
        end
        ccall(:jl_set_task_tid, Cvoid, (Any, Cint), task, tid - 1)
        schedule(task)
        push!(tasks, task)
    end
    foreach(wait, tasks)
    T = typeof(time_ns())
    notified = T[r[] for r in results if r[] isa Number]
    pre, post = only(r[] for r in results if !(r[] isa Number))
    return (; notified, pre, post)
end

function main(args = ARGS)
    output, = args
    mkpath(dirname(output))
    data = Dict(
        :results => sample(notify_overheads, 1000),
        :metadata => metadata(),
    )
    open(output, write = true) do io
        JSON.print(io, data)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
"""
Aircraft object containing state variables and parameters.
"""
Base.@kwdef mutable struct Aircraft
	r::Vector{R} = zeros(3) # position [ft]
	v::Vector{R} = zeros(3) # velocity [ft/s]

	# constants
	a_max::R     = G/2      # maximum acceleration along any axis [ft/s^2]
	v_min::R     = 50.0     # minimum velocity [ft/s]
	z_min::R     = 1.0      # minimum altitude [ft]
end

vec(ac::Aircraft) = [ac.r; ac.v]

"""
Command object specifying input to aircraft controls. Typically generated by transition
distribution.
"""
Base.@kwdef mutable struct Command
	h_d::R   = 0.0  # vertical rate [ft/s]
	vg_d::R  = 0.0  # ground acceleration [ft/s^2]
	psi_d::R = 0.0  # turn rate [deg/s]
end

Command(arr::Vector) = Command(arr...)
vec(cmd::Command) = [cmd.h_d; cmd.vg_d; cmd.psi_d]

function observation(ac::Aircraft; scale::Bool=true)
	r0 = scale ? 50.0 : 1.0
	v0 = scale ? 100.0 : 1.0
	obs = [ac.r / r0; ac.v / v0]
	return obs
end

function initialize(ac::Aircraft, ac_enc::Aircraft)
	ac.r = deepcopy(ac_enc.r)
	ac.v = deepcopy(ac_enc.v)
end

function update(ac::Aircraft, cmd::Command, dt)
    r, v = ac.r, ac.v
    vg = sqrt(v[1]^2 + v[2]^2)
	h_d, vg_d, psi_d = cmd.h_d, cmd.vg_d, deg2rad(cmd.psi_d)

	# minimum velocity
	vg_d = vg < ac.v_min ? max(vg_d, 0.0) : vg_d

	# compute accelerations
    ax = vg_d * v[1] / vg + v[2] * psi_d
    ay = vg_d * v[2] / vg - v[1] * psi_d
    az = (h_d - v[3]) / dt
    a = [ax; ay; az]
    clamp!(a, -ac.a_max, ac.a_max)

	# update state
    v_new = v + a * dt
    r_new = r + (v + v_new) * dt / 2

	# minimum altitude
	if r_new[3] < ac.z_min
		r_new[3] = ac.z_min
		v_new[3] = 0.0
	end

    ac.r, ac.v = r_new, v_new
end

distance(ac1::Aircraft, ac2::Aircraft) = norm(ac1.r - ac2.r)

"""
Return tuple of horizontal and vertial separation, optionally scaled by NMAC parameters.
"""
function separation(ac1::Aircraft, ac2::Aircraft; scale::Tuple{R,R}=(1.0, 1.0))
	r1, r2 = ac1.r, ac2.r
	return (norm(r1[1:2] - r2[1:2]), (scale[1] / scale[2]) * abs(r1[3] - r2[3]))
end

"""
Reduce symmetries of two-aircraft system, yielding R^12 -> R^8. Transformed coordinate space
is in a rotated frame of reference where the center of the system is at the origin. The
aircraft lie along the x-axis (except for displacement along the z-axis) with aircraft 1 at
position (d, 0, z) and aircraft 2 at (-d, 0, -z). Horizontal velocity is expressed as a
magnitude and angle is measured clockwise from the +y direction for both aircraft.

Relative positions are further transformed with an arcsinh function. This provides log-like
squashing at large magnitudes, preserves linear behavior at small magnitudes, and extends
smoothly and antisymmetrically to negative values. Since aircraft may travel large distances
before collision, the transformation encourages the learning of scale-dependent features.
"""
function kernel(ac1::Aircraft, ac2::Aircraft; scale::Bool=true)
	ac1, ac2 = observation.((ac1, ac2); scale=scale)
	r1, v1 = ac1[1:3], ac1[4:6]
	r2, v2 = ac2[1:3], ac2[4:6]

	rc = (r1 + r2) / 2        # center of two-aircraft system
	dr = r1 - rc              # ac1 position relative to center
	d = norm(dr[1:2])         # horizontal separation
	phi = atan(dr[2], dr[1])  # angle of two-aircraft system
	vg1 = norm(v1[1:2])       # ac1 ground speed
	vg2 = norm(v2[1:2])       # ac2 ground speed
	psi1 = atan(v1[1:2]...)   # ac1 heading
	psi2 = atan(v2[1:2]...)   # ac2 heading

	state = zeros(8)
	# 2D components
	state[1] = asinh(d)
	state[2] = wrap_to_pi(psi1 + phi) # ac1 relative heading
	state[3] = vg1
	state[4] = wrap_to_pi(psi2 + phi) # ac2 relative heading
	state[5] = vg2
	# 3D components
	state[6] = asinh(dr[3])
	state[7] = v1[3]
	state[8] = v2[3]
	return state
end


text:
hes_reshape(h, len) = [h[i,j] for i=1:len for j=1:i]

function hes_structure!(rows, cols, len)
  idx = 1
  for row in 1:len
    for col in 1:row
      rows[idx] = row
      cols[idx] = col
      idx += 1
    end
  end
  nothing
end

function ipopt_hessian(xarg, mode, rows, cols, obj_factor, lams, values,
                       hessfn, constr_hessv, nconstr)
  (mode == :Structure) && return hes_structure!(rows, cols, length(xarg))
  @assert length(lams) == nconstr "Disagreement in lengths of lambdas and constraint functions."
  h = hessfn(xarg)
  values .= hes_reshape(h, length(xarg)).*obj_factor
  for (lj, hesconstrj) in zip(lams, constr_hessv)
    constrj_hes = hesconstrj(xarg)
    values .+= lj*hes_reshape(constrj_hes, length(xarg))
  end
end

function jac_structure!(rows, cols, len, nconstr)
  for (j, ix) in enumerate(Iterators.partition(1:(len*nconstr), len))
    rows[ix].=j
    cols[ix].=collect(1:len)
  end
  nothing
end

function ipopt_constr_jac(xarg, mode, rows, cols, values, g_jac, nconstr)
  (mode == :Structure) && return jac_structure!(rows, cols, length(xarg), nconstr)
  values .= g_jac(xarg)
end


text:
# This file investigates integration error with respect to floating point precison.

using FractalTools
using Makie

# Construct interpolation data
f(x, y) = x^2 + y^2 + 1
precisions = 2 .^ (7 : 11)

# Specify theoretical value
fval = 35 / 48

npts = 100
mse = map(precisions) do prec
    # Set precision
    setprecision(prec) do
        freevar = BigFloat(0.001)

        # Construct interpolation data
        ngon = [
            BigFloat.([0.0, 0.0]),
            BigFloat.([1.0, 0.0]),
            BigFloat.([0.5, 1.0])
        ]
        pts = getdata(f, ngon, npts)

        # Compute integration.
        ival = integrate(pts, Interp2D(freevar))

        # Compute errro
        abs(fval - ival)
    end
end

# Plot mse
fig = Figure()
ax = fig[1, 1] = Axis(fig, xlabel="Precision", ylabel="Integration Error", title="2D Integration Error")
stem!(ax, precisions, mse, color=:black)
save(joinpath(@__DIR__, "integ2d_error.png"), fig)
display(fig)


text:
using Base.Iterators, LinearAlgebra, Optim, Plotly;
import DelimitedFiles.readdlm;

NDimensions = Int64(4);
#Define the single qubit measurements, R is the state with "-i".
H = Float64[1, 0];
V = Float64[0, 1];
statevars = [H, V];

Id = Matrix{Float64}(I, 2, 2);
Hx = [1 1; 1 -1]/sqrt(2);
Hy = [1 1; im -im]/sqrt(2);
Trafovars = [Id, Hx, Hy];

Trafos = [kron(l,k) for k in Trafovars, l in Trafovars];
StateArray = [kron(l,k) for k in statevars, l in statevars];
StateMeasured = [];

for a in 1:length(Trafos)
    for b in 1:NDimensions
        push!(StateMeasured,Trafos[a]*StateArray[b])
    end
end

D2Array = readdlm("D:\\myJulia\\two_qubit_tomo_measurement_2.txt", skipstart=2)[1:end,34:37];
Data = vec(transpose(D2Array));
SigmaData = sqrt.(Data);

function BadnessPolynomial!(t)
    pM = [t[1] 0 0 0;
        t[5] + im*t[6]  t[2] 0 0;
        t[11] + im*t[12] t[7] + im*t[8] t[3] 0;
        t[15] + im*t[16] t[13] + im*t[14] t[9] + im*t[10] t[4]];

    pMd = conj(transpose(pM));

    global GeneralDM = pMd*pM;
    Prediction(DM, State) = sum(conj(transpose(State))*DM*State);

    BD = sum((Prediction(GeneralDM, StateMeasured[i]) - Data[i])^2/Prediction(GeneralDM, StateMeasured[i]) for i=1:length(StateMeasured));
    return real(BD)
end

Bell = (kron(H,V)-kron(V,H))/sqrt(2);
IdealDM = reshape(kron(Bell, conj(Bell)),(NDimensions,NDimensions))

t0 = vec(IdealDM);
@time res = optimize(BadnessPolynomial!, t0, NelderMead(), Optim.Options(g_tol = 10e-10, iterations = 20000));
#@time res = optimize(BadnessPolynomial!, t0, LBFGS(), Optim.Options(iterations = 20000)) #Requires a function and gradient (will be approximated if omitted)
#inner_optimizer = GradientDescent()
#lower = fill(-1.0, NDimensions);
#upper = fill(1.0, NDimensions);
#@time res = optimize(BadnessPolynomial!, lower, upper, t0, Fminbox(inner_optimizer), Optim.Options(iterations = 50000)) #Box constrained minimization
temp = Optim.minimizer(res);

MaxLikDM = GeneralDM/tr(GeneralDM)

Fidelity = tr(sqrt(sqrt(MaxLikDM)*IdealDM*sqrt(MaxLikDM)))^2;
Purity = tr(MaxLikDM*MaxLikDM);
LinearEntropy = NDimensions/(NDimensions - 1)*(1 - Purity);

function plotme3d(dta)
    x = [];
    y = [];
    re_z = [];
    im_z = [];
    for i = 1:size(dta)[1]
        for j = 1:size(dta)[1]
            push!(x, i);
            push!(y, j);
            push!(re_z, real(dta[i,j]));
            push!(im_z, imag(dta[i,j]));
        end
    end

    trace1 = scatter3d(x = x, y = y, z = re_z, mode = "markers", marker = attr(size = 5), name = "Real part");
    trace2 = scatter3d(x = x, y = y, z = im_z, mode = "markers", marker = attr(size = 5), name = "Imaginary part");
    layout = Layout(xaxis = attr(tickvals=[1,2,3,4], ticktext = ["HH", "HV", "VH", "VV "]),
            yaxis = attr(tickvals=[1,2,3,4], ticktext = ["HH", "HV", "VH", "VV "]));
    #plot([trace1, trace2])
    p1 = Plotly.plot(trace1, layout);
    p2 = Plotly.plot(trace2, layout);
    [p1, p2]
end

plotme3d(MaxLikDM)


text:
# ---
# title: 1326. Minimum Number of Taps to Open to Water a Garden
# id: problem1326
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Greedy
# link: <https://leetcode.com/problems/minimum-number-of-taps-to-open-to-water-a-garden/description/>
# hidden: true
# ---
#
# There is a one-dimensional garden on the x-axis. The garden starts at the
# point `0` and ends at the point `n`. (i.e The length of the garden is `n`).
#
# There are `n + 1` taps located at points `[0, 1, ..., n]` in the garden.
#
# Given an integer `n` and an integer array `ranges` of length `n + 1` where
# `ranges[i]` (0-indexed) means the `i-th` tap can water the area `[i -
# ranges[i], i + ranges[i]]` if it was open.
#
# Return _the minimum number of taps_ that should be open to water the whole
# garden, If the garden cannot be watered return **-1**.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/01/16/1685_example_1.png)
#
#
#
#     Input: n = 5, ranges = [3,4,1,1,0,0]
#     Output: 1
#     Explanation: The tap at point 0 can cover the interval [-3,3]
#     The tap at point 1 can cover the interval [-3,5]
#     The tap at point 2 can cover the interval [1,3]
#     The tap at point 3 can cover the interval [2,4]
#     The tap at point 4 can cover the interval [4,4]
#     The tap at point 5 can cover the interval [5,5]
#     Opening Only the second tap will water the whole garden [0,5]
#
#
# **Example 2:**
#
#
#
#     Input: n = 3, ranges = [0,0,0,0]
#     Output: -1
#     Explanation: Even if you activate all the four taps you cannot water the whole garden.
#
#
# **Example 3:**
#
#
#
#     Input: n = 7, ranges = [1,2,1,0,2,1,0,1]
#     Output: 3
#
#
# **Example 4:**
#
#
#
#     Input: n = 8, ranges = [4,0,0,0,0,0,0,0,4]
#     Output: 2
#
#
# **Example 5:**
#
#
#
#     Input: n = 8, ranges = [4,0,0,0,4,0,0,0,4]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10^4`
#   * `ranges.length == n + 1`
#   * `0 <= ranges[i] <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end