text:
########################################################################
# Level set method
########################################################################
"""
Return an exit flag for the level-set oracle.
"""
function checkExitLvl(gap::Float64, ℓ::Float64, u::Float64, k::Int64, α::Float64, feaTol::Float64, maxIts::Int64, rule::String)
    if rule == "newton"
        # u < α + feaTol && return :feasible              # |Mx-b|^2/2 < α is nearly satisfied
        ℓ > α + feaTol && return :suboptimal            # ℓ is large enough to update τ
    elseif rule == "bisection"
        u < α && return :suboptimal_large               # τ too large
        ℓ > α && return :suboptimal_small               # τ too small
    end
    gap ≤ feaTol * eps(Float64) && return :optimal      # current lasso problem is fully optimal
    k ≥ maxIts && return :iterations                    # out of iterations
    return :noerror                                     # otherwise
end

"""
Dual conditional gradient method for the problem

    minimize 1/2||Mx-b||^2  subj to  gauge(x|A)≤τ

Includes various stopping conditions needed for
the level-set methods.
"""
function solveSubproblem!(dcg::DualCGIterable, τ::Float64, α::Float64, feaTol::Float64, maxIts::Int64, rule::String)
    setRadius!(dcg, τ)
    k = 0; u = 0.0; ℓ = 0.0
    exitFlag = :noerror
    for (gap, r) in dcg
        k = k + 1
        u = norm(r)^2/2
        ℓ = max(α, u - gap)
        exitFlag = checkExitLvl(gap, ℓ, u, k, α, feaTol, maxIts, rule)
        exitFlag == :noerror || break
    end
    s = support(dcg.A, dcg.z)
    return u, ℓ, s, k, exitFlag
end

"""
Level-set method for the problem

    min gauge(x | A) subj to 1/2||Mx-b||^2 ≤ α
"""
function level_set(M::AbstractLinearOp, b::Vector{Float64}, A::AbstractAtomicSet;
                   α::Float64 = 0.0,
                   tol::Float64 = 1e-12,
                   gapTol::Float64 = 1+norm(b),
                   maxIts::Int = size(M,2),
                   pr::Bool = true,
                   logger::Bool = true,
                   rule::String = "newton",
                   τmax::Float64 = 1.0,
                   dual_gaps::Vector{Float64} = Vector{Float64}(),
                   infeas::Vector{Float64} = Vector{Float64}()
                   )

    feaTol = (tol * norm(b))^2/2
    # Structure to hold the iterates
    dcg = DualCGIterable(M, b, A)
    # Structure to hold the solution
    sol = Solution(M, b, A)

    # head logging
    m, n = size(M)
    logger && logger_head_lvl(m, n, feaTol, α, maxIts)

    # Check whether b is zero
    if norm(b) < tol
        return sol
    end

    # Initialization at τ = 0
    τ = 0.0
    u = norm(b)^2/2; ℓ = norm(b)^2/2
    if rule == "newton"
        z = getNegGradient(dcg)
        s = support(A, z)
    elseif rule == "bisection"
        τmin = 0.0
    else
        println("invalid updating rule")
        return sol
    end

    # ----------------------------------------------------------------
    # min loop
    # ----------------------------------------------------------------
    k = 0; totIts = 0; exitFlag = :suboptimal_small
    while true
        k = k + 1

        # --------------------------------------------------------------
        # Update level-set radius.
        # --------------------------------------------------------------
        if rule == "newton"
            τ = τ + (ℓ-α)/s
        elseif rule == "bisection"
            if exitFlag == :suboptimal_small
                τmin = τ
            elseif exitFlag == :suboptimal_large
                τmax = τ
            end
            τ = (τmin + τmax)/2
            ρ = τmax - τmin
            @show ρ
        end

        # --------------------------------------------------------------
        # Check for exit.
        # --------------------------------------------------------------
        exitFlag == :iterations && break
        exitFlag == :feasible && break
        exitFlag == :optimal && break

        # --------------------------------------------------------------
        # Solve subproblem to obtain lower minorant defined by (ℓ,s).
        # --------------------------------------------------------------
        u, ℓ, s, minorIts, exitFlag = solveSubproblem!(dcg, τ, α, feaTol, maxIts-totIts, rule)

        # --------------------------------------------------------------
        # Primal recovery.
        # --------------------------------------------------------------
        # if pr && u - α ≤ gapTol
        if pr
            exitFlag = primalrecover!(dcg, sol, α, feaTol, exitFlag)
        end

        # --------------------------------------------------------------
        # record dual gaps and primal recover infeas
        # --------------------------------------------------------------
        push!(dual_gaps, dual_obj_gap(dcg, τmax, α))
        push!(infeas, sol.feasnew)

        # --------------------------------------------------------------
        # Logging and bookkeeping.
        # --------------------------------------------------------------
        logger && logger_level_lvl(α, τ, k, minorIts, ℓ, u, exitFlag, sol.feas)
        totIts += minorIts
    end

    primalrecover!(dcg, sol, α, feaTol, exitFlag)
    sort!(dual_gaps, rev=true)
    sort!(infeas, rev=true)

    # foot logging
    logger&& logger_foot_lvl(b, sol.feas, totIts, u-ℓ)

    return sol
end

function logger_head_lvl(m, n, feaTol, α, maxIts)
    println()
    println("  -------------------------------------------------------------------------")
    println("  Polar Level Set Method")
    println("  -------------------------------------------------------------------------")
    @printf("  %-22s %7d %7s %-20s %7d\n","number of variables",n,"","number of constraints",m)
    @printf("  %-22s %7.2e %7s %-20s %7.2e\n","feasibility tolerance",feaTol,"","α",α)
    @printf("  %-22s %7d \n","max iterations",maxIts)
    println("  -------------------------------------------------------------------------")
    @printf("  %5s   %8s   %8s   %8s   %8s   %8s  %15s  %10s\n",
            "Major","Minor","u-α","ℓ-α","gap","τ","infeas-α","Subproblem")
end

function logger_foot_lvl(b, feas, totIts, gap)
    scale = 1+norm(b)
    normr = sqrt(2*feas)
    println("  -------------------------------------------------------------------------")
    @printf("  %-22s %8.1e %7s %-22s %8.2e\n","residual (abs)",normr,"","optimality gap (abs)",gap)
    @printf("  %-22s %8.1e %7s %-22s %8.2e\n","residual (rel)",normr/scale,"","optimality gap (rel)",gap/scale)
    @printf("  %-22s %8d %7s %-22s %8s\n","total iterations",totIts,"","","")
    println("  -------------------------------------------------------------------------")
end

function logger_level_lvl(α, τ, k, minorItns, ℓ, u, exitFlag, feas)
    gap = u - ℓ
    # optimal = log10(dot(p.b,p.oracle.r)/norm(p.oracle.z,Inf))
    @printf("  %5d   %8d   %8.2e   %8.2e   %8.2e   %8.2e   %12.2e   %s\n",
            k, minorItns, u-α, ℓ-α, gap, τ, feas-α, exitFlag)
end

function dual_obj_gap(p::DualCGIterable, τ::Float64, α::Float64)
    β = τ
    s = support(p.A, p.z)
    y = deepcopy(p.r); y .*= 1 / s
    dobj = norm(y)^2/(2*β) + β*α - dot(y, p.b)
    return abs(dobj + τ)
end


text:
function Beta_tr(mu, sd)
    a = ((1 - mu) / sd ^ 2 - 1 / mu) * mu ^ 2
    b = a * (1 / mu - 1)
    return a, b
end

function Gamma_tr(mu, sd)
    b = sd ^ 2 / mu
    a = mu / b
    return a, b
end

function InvGamma_tr(mu, sd)
    a = mu ^ 2 / sd ^ 2 + 2
    b = mu * (a - 1)
    return a, b
end

function variance_check(u0)
    u0_variance = u0.C.U' * u0.C.U
    return (maximum(abs.(u0_variance)) > 1e10)
end

@model function rbc_kalman(z, m, p_f, α_prior, β_prior, ρ_prior, cache, settings)
    α ~ truncated(Normal(α_prior[1], α_prior[2]), α_prior[3], α_prior[4])
    β_draw ~ Gamma(β_prior[1], β_prior[2])
    ρ ~ Beta(ρ_prior[1], ρ_prior[2])
    β = 1 / (β_draw / 100 + 1)
    p_d = (; α, β, ρ)
    (settings.print_level > 0) && @show p_d
    T = size(z, 2)
    sol = generate_perturbation(m, p_d, p_f, Val(1); cache)
    (settings.print_level > 1) && println("Perturbation generated")

    if !(sol.retcode == :Success) || variance_check(sol.x_ergodic)
        (settings.print_level > 0) && println("Perturbation failed / Infinite variance with retcode $(sol.retcode)")
        @addlogprob! -Inf
    else
        (settings.print_level > 1) && println("Calculating likelihood")
        # Simulate and get the likelihood.
        problem = LinearStateSpaceProblem(sol.A, sol.B, sol.C, sol.x_ergodic, (0, T),
                                          noise = nothing, obs_noise = sol.D, observables = z)
        @addlogprob! solve(problem, KalmanFilter(); save_everystep = false).loglikelihood
    end
    return
end

@model function rbc_joint_1(z, m, p_f, α_prior, β_prior, ρ_prior, cache, settings, x0)
    α ~ truncated(Normal(α_prior[1], α_prior[2]), α_prior[3], α_prior[4])
    β_draw ~ Gamma(β_prior[1], β_prior[2])
    ρ ~ Beta(ρ_prior[1], ρ_prior[2])
    β = 1 / (β_draw / 100 + 1)
    p_d = (; α, β, ρ)
    (settings.print_level > 0) && @show p_d
    T = size(z, 2)
    ϵ_draw ~ MvNormal(m.n_ϵ * T, 1.0)
    ϵ = reshape(ϵ_draw, m.n_ϵ, T)
    sol = generate_perturbation(m, p_d, p_f, Val(1); cache)
    (settings.print_level > 1) && println("Perturbation generated")

    if !(sol.retcode == :Success)
        (settings.print_level > 0) && println("Perturbation failed / Infinite variance with retcode $(sol.retcode)")
        @addlogprob! -Inf
    else
        (settings.print_level > 1) && println("Calculating likelihood")
        # Simulate and get the likelihood.
        problem = LinearStateSpaceProblem(sol.A, sol.B, sol.C, x0, (0, T),
                                          noise = ϵ, obs_noise = sol.D, observables = z)
        @addlogprob! solve(problem, NoiseConditionalFilter(); save_everystep = false).loglikelihood
    end
    return
end

@model function rbc_joint_2(z, m, p_f, α_prior, β_prior, ρ_prior, cache, settings, x0)
    α ~ truncated(Normal(α_prior[1], α_prior[2]), α_prior[3], α_prior[4])
    β_draw ~ Gamma(β_prior[1], β_prior[2])
    ρ ~ Beta(ρ_prior[1], ρ_prior[2])
    β = 1 / (β_draw / 100 + 1)
    p_d = (; α, β, ρ)
    (settings.print_level > 0) && @show p_d
    T = size(z, 2)
    ϵ_draw ~ MvNormal(m.n_ϵ * T, 1.0)
    ϵ = reshape(ϵ_draw, m.n_ϵ, T)
    sol = generate_perturbation(m, p_d, p_f, Val(2); cache)
    (settings.print_level > 1) && println("Perturbation generated")

    if !(sol.retcode == :Success)
        (settings.print_level > 0) && println("Perturbation failed / Infinite variance with retcode $(sol.retcode)")
        @addlogprob! -Inf
    else
        (settings.print_level > 1) && println("Calculating likelihood")
        # Simulate and get the likelihood.
        problem = QuadraticStateSpaceProblem(sol.A_0, sol.A_1, sol.A_2, sol.B, sol.C_0, sol.C_1, sol.C_2, x0, (0, T),
                                             noise = ϵ, obs_noise = sol.D, observables = z)
        @addlogprob! solve(problem, NoiseConditionalFilter(); save_everystep = false).loglikelihood
    end
    return
end

@model function FVGQ20_kalman(z, m, p_f, params, cache, settings)
    T = size(z, 2)
    # Priors
    β_draw ~ Gamma(params.β[1], params.β[2])
    β = 1 / (β_draw / 100 + 1)
    h ~ Beta(params.h[1], params.h[2])
    κ ~ truncated(Normal(params.κ[1], params.κ[2]), params.κ[3], params.κ[4])
    χ ~ Beta(params.χ[1], params.χ[2])
    γR ~ Beta(params.γR[1], params.γR[2])
    γΠ ~ truncated(Normal(params.γΠ[1], params.γΠ[2]), params.γΠ[3], params.γΠ[4])
    Πbar_draw ~ Gamma(params.Πbar[1], params.Πbar[2])
    Πbar = Πbar_draw / 100 + 1
    ρd ~ Beta(params.ρd[1], params.ρd[2])
    ρφ ~ Beta(params.ρφ[1], params.ρφ[2])
    ρg ~ Beta(params.ρg[1], params.ρg[2])
    g_bar ~ Beta(params.g_bar[1], params.g_bar[2])
    σ_A ~ InverseGamma(params.σ_A[1], params.σ_A[2])
    σ_d ~ InverseGamma(params.σ_d[1], params.σ_d[2])
    σ_φ ~ InverseGamma(params.σ_φ[1], params.σ_φ[2])
    σ_μ ~ InverseGamma(params.σ_μ[1], params.σ_μ[2])
    σ_m ~ InverseGamma(params.σ_m[1], params.σ_m[2])
    σ_g ~ InverseGamma(params.σ_g[1], params.σ_g[2])
    Λμ ~ Gamma(params.Λμ[1], params.Λμ[2])
    ΛA ~ Gamma(params.ΛA[1], params.ΛA[2])
    # Likelihood
    θ = (; β, h, κ, χ, γR, γΠ, Πbar, ρd, ρφ, ρg, g_bar, σ_A, σ_d, σ_φ, σ_μ, σ_m, σ_g, Λμ, ΛA)
    (settings.print_level > 0) && @show θ
    sol = generate_perturbation(m, θ, p_f, Val(1); cache)
    (settings.print_level > 1) && println("Perturbation generated")

    if !(sol.retcode == :Success) || variance_check(sol.x_ergodic)
        (settings.print_level > 0) && println("Perturbation failed / Infinite variance with retcode $(sol.retcode)")
        @addlogprob! -Inf
    else
        z_trend = params.Hx * sol.x + params.Hy * sol.y
        z_detrended = z .- z_trend
        (settings.print_level > 1) && println("Calculating likelihood")

        # Simulate and get the likelihood.
        problem = LinearStateSpaceProblem(sol.A, sol.B, sol.C, sol.x_ergodic, (0, T),
                                          noise = nothing, obs_noise = sol.D, observables = z_detrended)
        @addlogprob! solve(problem, KalmanFilter(); save_everystep = false).loglikelihood
    end
    return
end

@model function FVGQ20_joint_1(z, m, p_f, params, cache, settings, x0)
    T = size(z, 2)
    # Priors
    β_draw ~ Gamma(params.β[1], params.β[2])
    β = 1 / (β_draw / 100 + 1)
    h ~ Beta(params.h[1], params.h[2])
    κ ~ truncated(Normal(params.κ[1], params.κ[2]), params.κ[3], params.κ[4])
    χ ~ Beta(params.χ[1], params.χ[2])
    γR ~ Beta(params.γR[1], params.γR[2])
    γΠ ~ truncated(Normal(params.γΠ[1], params.γΠ[2]), params.γΠ[3], params.γΠ[4])
    Πbar_draw ~ Gamma(params.Πbar[1], params.Πbar[2])
    Πbar = Πbar_draw / 100 + 1
    ρd ~ Beta(params.ρd[1], params.ρd[2])
    ρφ ~ Beta(params.ρφ[1], params.ρφ[2])
    ρg ~ Beta(params.ρg[1], params.ρg[2])
    g_bar ~ Beta(params.g_bar[1], params.g_bar[2])
    σ_A ~ InverseGamma(params.σ_A[1], params.σ_A[2])
    σ_d ~ InverseGamma(params.σ_d[1], params.σ_d[2])
    σ_φ ~ InverseGamma(params.σ_φ[1], params.σ_φ[2])
    σ_μ ~ InverseGamma(params.σ_μ[1], params.σ_μ[2])
    σ_m ~ InverseGamma(params.σ_m[1], params.σ_m[2])
    σ_g ~ InverseGamma(params.σ_g[1], params.σ_g[2])
    Λμ ~ Gamma(params.Λμ[1], params.Λμ[2])
    ΛA ~ Gamma(params.ΛA[1], params.ΛA[2])
    ϵ_draw ~ MvNormal(m.n_ϵ * T, 1.0)
    ϵ = reshape(ϵ_draw, m.n_ϵ, T)
    # Likelihood
    θ = (; β, h, κ, χ, γR, γΠ, Πbar, ρd, ρφ, ρg, g_bar, σ_A, σ_d, σ_φ, σ_μ, σ_m, σ_g, Λμ, ΛA)
    (settings.print_level > 0) && @show θ
    sol = generate_perturbation(m, θ, p_f, Val(1); cache)
    (settings.print_level > 1) && println("Perturbation generated")
    if !(sol.retcode == :Success)
        (settings.print_level > 0) && println("Perturbation failed / Infinite variance with retcode $(sol.retcode)")
        @addlogprob! -Inf
    else
        z_trend = params.Hx * sol.x + params.Hy * sol.y
        z_detrended = z .- z_trend
        # Simulate and get the likelihood.
        problem = LinearStateSpaceProblem(sol.A, sol.B, sol.C, x0, (0, T),
                                          noise = ϵ, obs_noise = sol.D, observables = z_detrended)
        @addlogprob! solve(problem, NoiseConditionalFilter(); save_everystep = false).loglikelihood
    end
    return
end

@model function FVGQ20_joint_2(z, m, p_f, params, cache, settings, x0)
    T = size(z, 2)
    # Priors
    β_draw ~ Gamma(params.β[1], params.β[2])
    β = 1 / (β_draw / 100 + 1)
    h ~ Beta(params.h[1], params.h[2])
    κ ~ truncated(Normal(params.κ[1], params.κ[2]), params.κ[3], params.κ[4])
    χ ~ Beta(params.χ[1], params.χ[2])
    γR ~ Beta(params.γR[1], params.γR[2])
    γΠ ~ truncated(Normal(params.γΠ[1], params.γΠ[2]), params.γΠ[3], params.γΠ[4])
    Πbar_draw ~ Gamma(params.Πbar[1], params.Πbar[2])
    Πbar = Πbar_draw / 100 + 1
    ρd ~ Beta(params.ρd[1], params.ρd[2])
    ρφ ~ Beta(params.ρφ[1], params.ρφ[2])
    ρg ~ Beta(params.ρg[1], params.ρg[2])
    g_bar ~ Beta(params.g_bar[1], params.g_bar[2])
    σ_A ~ InverseGamma(params.σ_A[1], params.σ_A[2])
    σ_d ~ InverseGamma(params.σ_d[1], params.σ_d[2])
    σ_φ ~ InverseGamma(params.σ_φ[1], params.σ_φ[2])
    σ_μ ~ InverseGamma(params.σ_μ[1], params.σ_μ[2])
    σ_m ~ InverseGamma(params.σ_m[1], params.σ_m[2])
    σ_g ~ InverseGamma(params.σ_g[1], params.σ_g[2])
    Λμ ~ Gamma(params.Λμ[1], params.Λμ[2])
    ΛA ~ Gamma(params.ΛA[1], params.ΛA[2])
    ϵ_draw ~ MvNormal(m.n_ϵ * T, 1.0)
    ϵ = reshape(ϵ_draw, m.n_ϵ, T)
    # Likelihood
    θ = (; β, h, κ, χ, γR, γΠ, Πbar, ρd, ρφ, ρg, g_bar, σ_A, σ_d, σ_φ, σ_μ, σ_m, σ_g, Λμ, ΛA)
    (settings.print_level > 0) && @show θ
    sol = generate_perturbation(m, θ, p_f, Val(2); cache)
    (settings.print_level > 1) && println("Perturbation generated")
    if !(sol.retcode == :Success)
        (settings.print_level > 0) && println("Perturbation failed / Infinite variance with retcode $(sol.retcode)")
        @addlogprob! -Inf
    else
        z_trend = params.Hx * sol.x + params.Hy * sol.y
        z_detrended = z .- z_trend
        # Simulate and get the likelihood.
        problem = QuadraticStateSpaceProblem(sol.A_0, sol.A_1, sol.A_2, sol.B, sol.C_0, sol.C_1, sol.C_2, x0, (0, T),
                                             noise = ϵ, obs_noise = sol.D, observables = z_detrended)
        @addlogprob! solve(problem, NoiseConditionalFilter(); save_everystep = false).loglikelihood
    end
    return
end

ext: jl

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples


text:
module Menus

Base.Experimental.@optlevel 2

using Redux
using CImGui
using ..MenuItems
import ..MenuItems: AbstractMenuItemAction, AbstractGenericMenuItem, AbstractMenuItem
using ..ToggleMenuItems
import ..ToggleMenuItems: ToggleMenuItem

## actions
"""
    AbstractGenericMenuAction <: AbstractSyncAction
Abstract supertype for all kinds of menu actions. Literally, its subtype could be anything.
"""
abstract type AbstractGenericMenuAction <: AbstractSyncAction end

"""
    AbstractMenuAction <: AbstractGenericMenuAction
Abstract supertype for all menu actions. Any subtypes of this type should has a
`label` field or at least `x.label` should work.
"""
abstract type AbstractMenuAction <: AbstractGenericMenuAction end

# for determining whether we should apply the reducer function to certain action types.
struct Valid end
struct Invalid end

is_valid(a::AbstractSyncAction) = Invalid
is_valid(a::AbstractGenericMenuAction) = Invalid
is_valid(a::AbstractMenuAction) = Valid

_get_label(a::AbstractSyncAction, ::Type{Invalid}) = "___REDUXIMGUI_RESERVED_DUMMY_ACTION_LABEL"
_get_label(a::AbstractSyncAction, ::Type{Valid}) = a.label
get_label(a::AbstractSyncAction) = _get_label(a, is_valid(a))

"""
    Rename(label, new_label)
Change widget's label to `new_label`.
Note that renaming may also change the identifier, please refer to `help?> CImGui.PushID`.
"""
struct Rename <: AbstractMenuAction
    label::String
    new_label::String
end

"""
    SetTriggeredTo(label, is_triggered)
Update widget's `is_triggered` state.
"""
struct SetTriggeredTo <: AbstractMenuAction
    label::String
    is_triggered::Bool
end

"""
    EditMenuItems(label, action::AbstractSyncAction)
Edit manu items.
"""
struct EditMenuItems{T<:AbstractSyncAction} <: AbstractMenuAction
    label::String
    action::T
end

"""
    Enable(label)
Enable the menu.
"""
struct Enable <: AbstractMenuAction
    label::String
end

"""
    Disable(label)
Disable the menu.
"""
struct Disable <: AbstractMenuAction
    label::String
end

"""
    AddMenu(label::AbstractString, items = [], is_enabled = true)
Add a [`Menu`](@ref).
"""
struct AddMenu <: AbstractMenuAction
    label::String
    items::Vector{AbstractImmutableState}
    is_enabled::Bool
end
AddMenu(label::AbstractString, items = []) = AddMenu(label, items, true)

"""
    DeleteMenu(label::AbstractString)
Delete the [`Menu`](@ref).
"""
struct DeleteMenu <: AbstractMenuAction
    label::String
end

## state
"""
    AbstractGenericMenu <: AbstractImmutableState
Abstract supertype for all kinds of menus. Literally, its subtype could be anything.
"""
abstract type AbstractGenericMenu <: AbstractImmutableState end

"""
    AbstractMenu <: AbstractGenericMenu
Abstract supertype for all menus. Any subtypes of this type should has a
`label` field or at least `x.label` should work.
"""
abstract type AbstractMenu <: AbstractGenericMenu end

is_valid(s::AbstractImmutableState) = Invalid
is_valid(s::AbstractGenericMenu) = Invalid
is_valid(s::AbstractMenu) = Valid

_get_label(s::AbstractImmutableState, ::Type{Invalid}) = "___REDUXIMGUI_RESERVED_DUMMY_STATE_LABEL"
_get_label(s::AbstractImmutableState, ::Type{Valid}) = s.label
get_label(s::AbstractImmutableState) = _get_label(s, is_valid(s))

"""
    Menu(label::AbstractString, items = [], is_enabled = true)
A menu state which contains a label a.k.a the identifier, a list of menu items,
a flag value `is_enabled` and a flag value `is_triggered` that records the state of
the latest poll events.
"""
struct Menu <: AbstractMenu
    label::String
    items::Vector{AbstractImmutableState}
    is_enabled::Bool
    is_triggered::Bool
end
Menu(label::AbstractString, items = [], is_enabled = true) =
    Menu(label, items, is_enabled, false)

is_enabled(s::Menu) = s.is_enabled
is_triggered(s::Menu) = s.is_triggered

# reducers
reducer(state::AbstractState, action::AbstractAction) = state
reducer(state::Vector{<:AbstractState}, action::AbstractAction) = state
reducer(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

reducer(s::Menu, a::Rename) = Menu(a.new_label, s.items, s.is_enabled, s.is_triggered)
reducer(s::Menu, a::Enable) = Menu(s.label, s.items, true, s.is_triggered)
reducer(s::Menu, a::Disable) = Menu(s.label, s.items, false, s.is_triggered)
reducer(s::Menu, a::SetTriggeredTo) = Menu(s.label, s.items, s.is_enabled, a.is_triggered)

function reducer(s::Menu, a::EditMenuItems)
    new_items = MenuItems.reducer(s.items, a.action)
    new_items = ToggleMenuItems.reducer(new_items, a.action)
    new_items = map(new_items) do item
        item isa Menu ? Menus.reducer(item, a.action) : item
    end
    Menu(s.label, new_items, s.is_enabled, s.is_triggered)
end

reducer(s::Dict{String,<:AbstractImmutableState}, a::AbstractMenuAction) =
    Dict(k => (get_label(v) == get_label(a) ? reducer(v, a) : v) for (k, v) in s)

reducer(s::Vector{<:AbstractImmutableState}, a::AbstractMenuAction) = map(s) do s
    get_label(s) == get_label(a) ? reducer(s, a) : s
end
reducer(s::Vector{<:AbstractImmutableState}, a::AddMenu) =
    [s..., Menu(a.label, a.items, a.is_enabled, false)]
reducer(s::Vector{<:AbstractImmutableState}, a::DeleteMenu) =
    filter(s -> s.label !== get_label(a), s)

# UI
"""
    (::Menu)(store::AbstractStore, get_state=Redux.get_state, chain_action=identity) -> Bool
Return `true` when triggered/activated.
`get_state` is a router function that tells how to find the target state from `store`.
`chain_action` is for chaining upstream actions.

To render anything as an MenuItem, you could use functors, for example,
```
struct ItemSeparator <: AbstractGenericMenuItem end
(x::ItemSeparator)() = CImGui.Separator()
```
"""
function (::Menu)(store::AbstractStore, get_state=Redux.get_state, chain_action=identity)
    s = get_state(store)
    is_activated = CImGui.BeginMenu(s.label, s.is_enabled)
    dispatch!(
        store,
        SetTriggeredTo(s.label, is_activated) |> chain_action
    )
    !is_activated && return false
    for (i, item) in enumerate(s.items)
        if item isa Union{AbstractMenuItem, AbstractMenu}
            item(
                store,
                x->get_state(x).items[i],
                x->EditMenuItems(s.label, x) |> chain_action,
            )
        elseif item isa Union{AbstractGenericMenuItem, AbstractGenericMenu}
            item()
        else
            CImGui.Text("menu item type not supported! skip rendering...")
        end
    end
    CImGui.EndMenu()
    return true
end

end # module


text:
"""
    mps!(u, domain, a, b, N; num_boundary = 2N, num_interior = 2N, optim_prec = 20)

Return the `λ` minimizing `σ(λ)` on the interval `[a, b]` and set `u`
to the corresponding eigenfunction.

The number of coefficients used for `u` is determined by `N`. The
number of boundary points and interior points used in the calculations
are given by `num_boundary` and `num_interior` respectively.

The minimum is computer to `optim_prec` bits of precision.

If `show_progress` is true then show a progress bar for the
minimization. If `show_trace` is true then show the trace of the
minimization. If `store_trace` is true also return the result of the
search for the minimum, if `a == b` (in which case no search is
performed) return `nothing`.
"""
function mps!(
    u::AbstractEigenfunction,
    domain::AbstractDomain,
    a::Real,
    b::Real,
    N::Integer = 8;
    num_boundary::Integer = 2N,
    num_interior::Integer = 2N,
    optim_prec::Integer = 20,
    show_progress::Bool = false,
    show_trace::Bool = false,
    store_trace::Bool = false,
)
    # Some types of eigenfunctions (e.g.
    # StandaloneLightningEigenfunction) needs that the number of terms
    # used is determined from the beginning.
    set_eigenfunction!(u, zeros(N))

    a > b && throw(ArgumentError("must have a < b, got a = $a, b = $b"))

    if a != b
        # Compute minimum of σ(λ)
        σ(λ) = sigma(λ, domain, u, N; num_boundary, num_interior)

        tol = 2.0^(-(optim_prec - 1))

        if show_progress
            start_prec = -log2((b - a) / ((b + a) / 2))

            optim_progress(state) = begin
                if state isa AbstractArray
                    state = state[end]
                end
                if state == :done
                    progress = "done"
                else
                    abs_error = state.metadata["x_upper"] - state.metadata["x_lower"]
                    rel_prec = -log2(abs_error / state.metadata["minimizer"])
                    progress = Float64(
                        max((rel_prec - start_prec) / (optim_prec - start_prec), 0),
                    )
                end
                @info "Computing minimum of σ(λ)" progress = progress
                false
            end
        else
            optim_progress = nothing
        end

        res = optimize(
            σ,
            a,
            b,
            rel_tol = tol,
            abs_tol = tol,
            callback = optim_progress;
            show_trace,
            store_trace,
        )

        if show_progress
            optim_progress(:done)
        end

        if !Optim.converged(res)
            @warn "Failed to compute minimum of σ(λ)"
        end

        λ = res.minimizer
    else
        λ = a
    end

    # Compute the eigenfunction corresponding to the computed minimum
    # of σ(λ).
    coefficients = sigma_coefficients(λ, domain, u, N; num_boundary, num_interior)
    set_eigenfunction!(u, coefficients)

    if store_trace
        if a != b
            return λ, res
        else
            return λ, nothing
        end
    else
        return λ
    end
end

"""
    iteratemps(domain, u, enclosure, Ns; kwargs...)

For each `N` in `Ns` compute the minimizing `λ` in `enclosure` and
corresponding eigenfunction with [`mps`](@ref). At each step also
compute an enclosure for `λ` and if this is better than the original
enclosure use that in the next step instead.

It takes a number of key-word arguments for handling progressively
updating parameters.
- `num_boundary_factor`, `num_interior_factor`: The number of points
  used on the boundary respectively in the interior is given by these
  factors multiplied by `N`.

There are a number of different ways for determining the tolerance to
use in the optimization of `σ(λ)` in the next step. It's
1. Constant: Use the tolerance every iteration. Set
`optim_prec_final` to the tolerance to use.
2. Linear: Specify a tolerance to use in the last iteration and then
linearly interpolate for the current `N`. Set `optim_prec_final` to
the tolerance and `optim_prec_linear` to true.
3. Step: Compute the current precision of the enclosure and use that
plus a fixed number as the tolerance for the next iteration. Set
`optim_prec_step` to true.
4. Adaptive: Estimate the improvement for the last step by computing
the precision of the current enclosure and the one in the iteration
before. Add the improvement to the precision of the current enclosure
plus a small constant. Set `optim_prec_adaptive` to true and
`optim_prec_adaptive_extra` to the extra precision to use.

The precision used in the computations is determined by the tolerance
for the optimization plus `extra_prec`.

If `recompute` is true then run `recompute!(u)` for each iteration.
This is useful for eigenfunction which precompute some of there
values, they can then be recomputed with the new precision. For most
eigenfunctions this doesn't have an effect.

If `rigorous_enclosure` is true then rigorously compute the enclosure
of the eigenvalue by rigorously bounding the maximum on the boundary
and the norm. If `rigorous_enclosure` is true then `rigorous_norm` can
be set to false to only compute the maximum rigorously.
"""
function iteratemps(
    domain::AbstractDomain,
    u::AbstractEigenfunction,
    enclosure::arb,
    Ns;
    num_boundary_factor::Integer = 2,
    num_interior_factor::Integer = 2,
    optim_prec_final = nothing,
    optim_prec_linear = false,
    optim_prec_step = nothing,
    optim_prec_adaptive = false,
    optim_prec_adaptive_extra = 10,
    extra_prec = 20,
    recompute = true,
    rigorous_enclosure = true,
    rigorous_norm = rigorous_enclosure,
    show_progress = false,
    show_trace = false,
    store_trace = false,
    extended_trace = false,
)
    # Make sure that at least one method for computing the precision
    # to use for the minimization is specified
    if isnothing(optim_prec_final) && isnothing(optim_prec_step) && !optim_prec_adaptive
        throw(
            ArgumentError(
                "at least one method for determining precision must be specified",
            ),
        )
    end

    # Setup function used for computing precision to which σ(λ) should
    # be minimized
    get_optim_prec(N, enclosure, λs, i) = begin
        optim_prec = typemax(Int)

        # If a final precision for the minimization is specified use
        # that, possibly with a linear interpolation.
        if !isnothing(optim_prec_final)
            if optim_prec_linear
                optim_prec = min(optim_prec, ceil(Int, optim_prec_final * N / Ns[end]))
            else
                optim_prec = min(optim_prec, optim_prec_final)
            end
        end

        # Set the precision to the relative accuracy of the current
        # enclosure plus a fixed precision
        if !isnothing(optim_prec_step)
            optim_prec =
                min(optim_prec, ArbTools.rel_accuracy_bits(enclosure) + optim_prec_step)
        end

        # Compute the precision increase between the last two
        # iterations and use the same plus a small constant
        if optim_prec_adaptive
            if i > 3 && isfinite(λs[i-2]) && isfinite(λs[i-1])
                optim_prec = min(
                    optim_prec,
                    max(
                        2ArbTools.rel_accuracy_bits(λs[i-1]) -
                        ArbTools.rel_accuracy_bits(λs[i-2]),
                        0,
                    ) + optim_prec_adaptive_extra,
                )
            else
                optim_prec = min(optim_prec, precision(domain.parent) - extra_prec)
            end
        end

        return optim_prec
    end

    # Copy the domain and eigenfunction to avoid changing the originals
    domain = deepcopy(domain)

    λs = similar(Ns, arb)
    us = [deepcopy(u) for _ in eachindex(Ns)]

    trace = MPSTrace()
    tracing = store_trace || show_trace || extended_trace

    if show_trace
        show(trace)
    end

    try
        for (i, N) in enumerate(Ns)
            ### Number of boundary and interior points to use ###
            num_boundary = num_boundary_factor * N
            num_interior = num_interior_factor * N

            ### Precision for optimization ###
            optim_prec = get_optim_prec(N, enclosure, λs, i)

            ### Precision to use in the computations ###
            new_prec = max(53, optim_prec + extra_prec) # Always use at least 53 bits

            ### Update the domain and eigenfunction with the new precision ###
            RR = RealField(new_prec)
            domain = typeof(domain)(domain, RR)
            set_domain!(us[i], domain)
            if recompute
                recompute!(us[i])
            end

            ### Run the computations ###
            @timeit_debug "$N" begin
                @timeit_debug "mps" begin
                    λ = setprecision(BigFloat, new_prec) do
                        mps!(
                            us[i],
                            domain,
                            getinterval(BigFloat, enclosure)...,
                            N,
                            store_trace = extended_trace;
                            num_boundary,
                            num_interior,
                            optim_prec,
                            show_progress,
                        )
                    end

                    if extended_trace
                        λ, optim_res = λ
                    end

                    λ = RR(λ)
                end

                @timeit_debug "enclosure" begin
                    if rigorous_enclosure
                        new_enclosure = enclose_eigenvalue(
                            domain,
                            us[i],
                            λ,
                            store_trace = tracing;
                            rigorous_norm,
                            extended_trace,
                            show_progress,
                        )
                    else
                        new_enclosure = enclose_eigenvalue_approx(
                            domain,
                            us[i],
                            λ,
                            max_numpoints = 4num_boundary_factor * N,
                            norm_numpoints = 4num_interior_factor * N,
                            store_trace = tracing;
                            extended_trace,
                        )
                    end
                end
            end

            if tracing
                metadata = Dict()

                if extended_trace
                    new_enclosure, n, m, maximize_trace = new_enclosure
                    metadata["optim_res"] = optim_res
                    metadata["maximize_trace"] = maximize_trace
                else
                    new_enclosure, n, m = new_enclosure
                end

                state = MPSState(N, new_prec, optim_prec, n, m, λ, new_enclosure, metadata)

                update!(trace, state, store_trace, show_trace)
            end

            λs[i] = new_enclosure
            if isfinite(new_enclosure)
                enclosure = setintersection(enclosure, new_enclosure)
            end
        end
    catch e
        if !(e isa InterruptException)
            rethrow(e)
        end
    end

    if store_trace
        return λs, us, trace
    else
        return λs, us
    end
end


text:
function yaml_node_from_string(fn::AbstractString)

    # Didn't find how to access the top node more easily.
    #
    yml_types = Dict{String,Function}()
    yml_types["!Cartesian"] = ((c,n)->n)
    yml_types["!Smolyak"] = ((c,n)->n)
    yml_types["!Normal"] = ((c,n)->n)
    yml_types["!MarkovChain"] = ((c,n)->n)
    yml_types["!Product"] = ((c,n)->n)
    yml_types["!PoissonProcess"] = ((c,n)->n)
    yml_types["!DeathProcess"] = ((c,n)->n)
    yml_types["!AgingProcess"] = ((c,n)->n)
    yml_types["!VAR1"] = ((c,n)->n)
    yml_types["tag:yaml.org,2002:null"] = ((c,n)->n)
    yml_types["tag:yaml.org,2002:bool"] = ((c,n)->n)
    yml_types["tag:yaml.org,2002:int"] = ((c,n)->n)
    yml_types["tag:yaml.org,2002:float"] = ((c,n)->n)
    yml_types["tag:yaml.org,2002:binary"] = ((c,n)->n)
    yml_types["tag:yaml.org,2002:timestamp"] = ((c,n)->n)
    yml_types["tag:yaml.org,2002:omap"] = ((c,n)->n)
    yml_types["tag:yaml.org,2002:pairs"] = ((c,n)->n)
    yml_types["tag:yaml.org,2002:set"] = ((c,n)->n)
    yml_types["tag:yaml.org,2002:str"] = ((c,n)->n)
    yml_types["tag:yaml.org,2002:seq"] = ((c,n)->n)
    yml_types["tag:yaml.org,2002:map"] = ((c,n)->n)

    return YAML.load(fn, yml_types)
end

function yaml_node_from_file(fn::AbstractString)
    txt = open(f->read(f,String), fn)
    txt = replace(txt, "\r"=>"")
    return yaml_node_from_string(txt)
end

mutable struct Location
    start_mark::YAML.Mark
    end_mark::YAML.Mark
end

mutable struct MissingElement <: Exception
    path::Vector{AbstractString}
    k::Integer # index of first missing element
    loc::Location
end

function Base.getindex(d::YAML.MappingNode, l::AbstractString...)
    el = d
    for k=1:length(l)
        # if !(typeof(el)<:YAML.MappingNode)
        #     println("This is not going to fly")
        if !(l[k] in keys(el))
            loc = Location(el.start_mark, el.end_mark)
            throw(MissingElement(collect(l), k, loc))
        else
            el = el[l[k]]
        end
    end
    return el
end
function Base.getindex(d::YAML.MappingNode, l::Symbol...)
    return Base.getindex(d, [string(e) for e in l]...)
end

mutable struct LinterException <: Exception
    msg::AbstractString
    loc::Location
    src::AbstractString
end

mutable struct LinterWarning <: Exception
    errvalue::AbstractString
    errtype::AbstractString
    msg::AbstractString
    loc::Location
    src::AbstractString
end

#LinterException(msg::AbstractString, loc::Location) = LinterException(msg, loc, "<string>")

LinterWarning(msg::AbstractString, loc::Location) = LinterWarning(msg, loc, "<string>", errtype)

function get_loc(d)
    Location(d.start_mark, d.end_mark)
end

function check_name(d)
  if !("name" in keys(d)) || typeof(d["name"].value ) != String || d["name"].value == ""
    return false
  end
end

function check_symbol_validity(sym)
  # Not sure if it is necesarry to check if that's a string: yaml structure saves by default as string(??)
  if  typeof(sym) != String || match(r"^[a-zA-Z0-9_]*$",sym) == nothing || string(sym)[1:1] == "_"
    return false
  end
end

#function check_equations(d::YAML.MappingNode, filename="<string>")
#Known_equation_types -> ["transition", "arbitrage", "value", "felicity", "expectation"]
#end

function check_model_sections(d::YAML.MappingNode, filename="<string>")
  errors = LinterWarning[]
  warnings = LinterWarning[]

  # Can add other reauired models sections: symbols, equations ...
  if check_name(d) == false
    if !("name" in keys(d))
      msg = ""
      errvalue = "name"
      errtype = "Missing model part"
      loc = get_loc(d)
    elseif typeof(d["name"].value ) != String || d["name"].value == ""
      errvalue = string(d["name"].value)
      errtype = "Invalid model name"
      loc = get_loc(d["name"])
      msg = "Model name should be non-empty string"
          # push!(errors/warnings, LinterWarning(errvalue, errtype, msg, loc, src)
    end

    push!(errors, LinterWarning(errvalue, errtype, msg, loc, filename))

  end

  return [errors, warnings]
end

function check_calibration(d::YAML.MappingNode, filename="<string>")
  errors = LinterWarning[]
  warnings = LinterWarning[]

  model_symbol_types = keys(d[:symbols])
  model_symbols = []
  model_symbols_vec = []
  for key in keys(d[:symbols])
      n = length(d[:symbols][key].value)
      for i=1:n
          # would be cooler to do directly for sym in syms[vg]
          sym = d[:symbols][key].value[i]
          push!(model_symbols_vec, (key, sym.value, sym))
          push!(model_symbols, sym.value)
      end
  end
  model_symbols = cat(model_symbols, keys(d[:definitions]); dims=1)

  for (i, key) in enumerate(keys(d["calibration"]))
    if !(key in model_symbols)

      errvalue = key
      errtype = "Calibration not in symbols"
      loc = get_loc(d["calibration"].value[i][1])
      msg = string(key, " is not declared in the symbols section ")
      push!(warnings, LinterWarning(errvalue, errtype, msg, loc, filename))

    elseif count(c -> c == key, collect( keys(d["calibration"]))) > 1
      floc = something(findfirst(isequal(key), keys(d["calibration"])), 0)
      loc_first = get_loc(d["calibration"].value[floc][1])
      if i > floc
        errvalue = key
        errtype = "Invalid calibration"
        loc = get_loc(d["calibration"].value[i][1])
        msg = string(key, " already declared in calibration section at line ",loc_first.start_mark.line, ", column ",loc_first.start_mark.column)
        push!(errors, LinterWarning(errvalue, errtype, msg, loc, filename))

      end
    end

  end

  return [errors, warnings]
end

function check_equations(d::YAML.MappingNode, filename="<string>")
  errors = LinterWarning[]
  warnings = LinterWarning[]

  # check symbol names:
  required_equation_types = ["transition"]
  optional_equation_types = ["arbitrage", "value", "felicity", "expectation", "direct_response"]
  known_equation_types = cat(required_equation_types, optional_equation_types; dims=1)

  model_equation_types = keys(d[:equations])

  node_transition = YAML.ScalarNode[]
  node_arbitrage = YAML.ScalarNode[]
  for eq in d[:equations].value
    if eq[1].value == "transition"
      node_transition = eq[1]
    elseif eq[1].value == "arbitrage"
      node_arbitrage = eq[1]
    end

  end

  for s in required_equation_types
      if !(s in model_equation_types)
          errvalue = string(s)
          errtype = "Missing equation type"
          msg = ""
          loc = get_loc(d[:equations])
          # push!(errors/warnings, LinterWarning(errvalue, errtype, msg, loc, src)
          push!(errors, LinterWarning(errvalue, errtype, msg, loc, filename))
      end
  end

  for (i,sg) in enumerate(keys(d["equations"]))
      if !(sg in cat(known_equation_types; dims=1))
          errvalue = string(sg)
          errtype = "Unknown equation type"
          msg = ""
          # get precise location
          loc = get_loc(d["equations"].value[i][1])
          push!(warnings, LinterWarning(errvalue, errtype, msg, loc, filename))
      end
  end

  if ("transition" in model_equation_types) && ("states" in keys(d[:symbols]))
    n_transition = length(d[:equations]["transition"].value)
    n_states = length(d[:symbols]["states"].value)

    if n_transition != n_states

      errvalue = string(node_transition.value)
      errtype = "Invalid number of equations"
      msg = "Number of transition equations should be equal to number of states"
      # get precise location
      loc = get_loc(node_transition)

      push!(errors, LinterWarning(errvalue, errtype, msg, loc, filename))
    end

  end

  if ("arbitrage" in model_equation_types) && ("controls" in keys(d[:symbols]))
    n_arbitrage = length(d[:equations]["arbitrage"].value)
    n_controls = length(d[:symbols]["controls"].value)

    if n_arbitrage != n_controls
      errvalue = string(node_arbitrage.value)
      errtype = "Invalid number of equations"
      msg = "Number of arbitrage equations should be equal to number of controls"
      # get precise location
      loc = get_loc(node_arbitrage.value)
      push!(warnings, LinterWarning(errvalue, errtype, msg, loc, filename))
    end

  end

  return [errors, warnings]
end

function check_symbols(d::YAML.MappingNode, filename="<string>")

    ### so far, this is just an example

    errors = LinterWarning[]
    warnings = LinterWarning[]

    # check symbol names:
    required_symbol_types = ["states", "controls", "exogenous", "parameters"]
    optional_symbol_types = [ "values", "rewards", "expectations"]
    known_symbol_types = cat(required_symbol_types, optional_symbol_types; dims=1)

    model_symbol_types = keys(d[:symbols])
    model_symbols = []
    model_symbols_vec = []
    for key in keys(d[:symbols])
        n = length(d[:symbols][key].value)
        for i=1:n
            # would be cooler to do directly for sym in syms[vg]
            sym = d[:symbols][key].value[i]
            push!(model_symbols_vec, (key, sym.value, sym))
            push!(model_symbols, sym.value)
        end
    end

    for s in required_symbol_types
        if !(s in model_symbol_types)
            errvalue = string(s)
            errtype = "Missing symbol type"
            msg = ""
            loc = get_loc(d[:symbols])
            # push!(errors/warnings, LinterWarning(errvalue, errtype, msg, loc, src)
            push!(errors, LinterWarning(errvalue, errtype, msg, loc, filename))
        end
    end

    for (i,sg) in enumerate(keys(d["symbols"]))
        if !(sg in cat(known_symbol_types; dims=1))
            errvalue = string(sg)
            errtype = "Unknown symbol type"
            msg = ""
            # get precise location
            loc = get_loc(d["symbols"].value[i][1])
            push!(warnings, LinterWarning(errvalue, errtype, msg, loc, filename))
        end
    end

    for (i,m) in enumerate(values(d["symbols"]))
      for (j, n) in enumerate(values(d["symbols"])[i])
        sym =n.value
        if check_symbol_validity(sym) == false
          errvalue = string(sym)
          errtype = "Invalid symbol"
          loc = get_loc(n)
          # get precise location
          if  typeof(sym) != String
            msg = "symbol should be a string"
          elseif match(r"^[a-zA-Z0-9_]*$",sym) == nothing
            msg = "symbol should be an 'alphanumeric' string"
          elseif string(sym)[1:1] == "_"
            msg = "symbol should not start with a number or an underscore"
          end
          push!(errors, LinterWarning(errvalue, errtype, msg, loc, filename))
        end
      end
    end

    for (i,symnode) in enumerate(model_symbols_vec)
      if count(c -> c == symnode[2], collect(model_symbols)) > 1
        floc = something(findfirst(isequal(symnode[2]), model_symbols), 0)
        loc_first = get_loc(model_symbols_vec[floc][3])
        if i > floc
          errvalue = symnode[2]
          errtype = "Invalid symbol"
          loc = get_loc(symnode[3])
          msg = string(symnode[2], " already declared as '", symnode[1] , " at line ",loc_first.start_mark.line, ", column ",loc_first.start_mark.column)
          push!(errors, LinterWarning(errvalue, errtype, msg, loc, filename))
        end
      end
    end

    for (i, sym) in enumerate(keys(d["definitions"]))

      if check_symbol_validity(sym) == false
        errvalue = string(sym)
        errtype = "Invalid definition"
        loc = get_loc(d["definitions"].value[i][1])
        # get precise location
        if  typeof(sym) != String
          msg = "symbol should be a string"
        elseif match(r"^[a-zA-Z0-9_]*$",sym) == nothing
          msg = "symbol should be an 'alphanumeric' string"
        elseif string(sym)[1:1] == "_"
          msg = "symbol should not start with a number or an underscore"
        end
        push!(errors, LinterWarning(errvalue, errtype, msg, loc, filename))
      end

    end

    for (i, key) in enumerate(keys(d["definitions"]))
      if (key in model_symbols)
        floc = something(findfirst(isequal(key), model_symbols), 0)
        loc_first = get_loc(model_symbols_vec[floc][3])
        declared_sym_type = model_symbols_vec[floc][1]

        errvalue = key
        errtype = "Invalid definition"
        loc = get_loc(d["definitions"].value[i][1])
        msg = string(key, "  already declared in symbols section as '", declared_sym_type , "' at line",loc_first.start_mark.line, ", column ",loc_first.start_mark.column)
        push!(errors, LinterWarning(errvalue, errtype, msg, loc, filename))

      elseif count(c -> c == key, collect( keys(d["definitions"]))) > 1
        floc = something(findfirst(isequal(key), keys(d["definitions"])), 0)
        loc_first = get_loc(d["definitions"].value[floc][1])
        if i > floc
          errvalue = key
          errtype = "Invalid definition"
          loc = get_loc(d["definitions"].value[i][1])
          msg = string(key, " already declared in definitions section at line ",loc_first.start_mark.line, ", column ",loc_first.start_mark.column)
          push!(errors, LinterWarning(errvalue, errtype, msg, loc, filename))

        end
      end
    end

    for (i,symnode) in enumerate(model_symbols_vec)
      if !(symnode[2] in keys(d["calibration"]))
          errvalue = symnode[2]
          errtype = "Symbol not calibrated"
          loc = get_loc(symnode[3])
          msg = string(symnode[2], " not found in calibration section")
          push!(warnings, LinterWarning(errvalue, errtype, msg, loc, filename))
      end
    end

    return [errors, warnings]

end

function check_equations(fn::AbstractString)
    d = yaml_node_from_file(fn)
    errs, wars = check_equations(d, fn)
end

function check_calibration(fn::AbstractString)
    d = yaml_node_from_file(fn)
    errs, wars = check_calibration(d, fn)
end

function check_symbols(fn::AbstractString)
    d = yaml_node_from_file(fn)
    errs, wars = check_symbols(d, fn)
end

function check_model_sections(fn::AbstractString)
    d = yaml_node_from_file(fn)
    errs, wars = check_model_sections(d, fn)
end

function check_model(fn::AbstractString)
    d = yaml_node_from_file(fn)
    errs, wars = check_model(d, fn)
end

function print_error(err::LinterWarning)
printstyled("error "; color=:light_red)
print("at line ",err.loc.start_mark.line, ", column ",err.loc.start_mark.column, " : ")
print(err.errtype)
if err.msg != ""
  printstyled(" '",err.errvalue,"' "; color=:light_green)
  println(">> ", err.msg)
else
  printstyled(" '",err.errvalue,"' "; color=:light_green)
  println()
end

end

function print_warning(err::LinterWarning)
printstyled("warning "; color=:light_blue)
print("at line ",err.loc.start_mark.line, ", colummn ",err.loc.start_mark.column, " : ")
print(err.errtype)
if err.msg != ""
  printstyled(" '",err.errvalue,"' "; color=:light_green)
  println(">> ", err.msg)
else
  printstyled(" '",err.errvalue,"' "; color=:light_green)
  println()
end
end

function format_human(errors::Vector{LinterWarning}, warnings::Vector{LinterWarning})
for err in cat(errors; dims=1)
    print_error(err)
end

for err in cat(warnings; dims=1)
    print_warning(err)
end
end

function lint(filename::AbstractString; format=:human)

  funnames = [check_model_sections, check_symbols, check_equations, check_calibration]
  errors = LinterWarning[]
  warnings = LinterWarning[]

  for fun in funnames
    errors2, warnings2 = fun(filename)
    errors = cat(errors, errors2; dims=1)
    warnings = cat(warnings, warnings2; dims=1)
  end

  if format == :human
      format_human(errors,warnings)
  else
      println("Format '", format, "' not implemented (yet).")
  end
  return cat(errors, warnings; dims=1)
end


text:
mutable struct MongoCollection
    _wrap_::Ptr{Cvoid}
    client::MongoClient
    db::AbstractString
    name::AbstractString

    MongoCollection(client::MongoClient, db::AbstractString, name::AbstractString) = begin
        dbCStr = String(db)
        nameCStr = String(name)
        collection = new(
            ccall(
                (:mongoc_client_get_collection, libmongoc),
                Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{UInt8}, Ptr{UInt8}),
                client._wrap_, dbCStr, nameCStr
                ),
            client,
            db,
            name
            )
        finalizer(collection, destroy)
        return collection
    end
end
export MongoCollection

function Base.show(io::IO, collection::MongoCollection)
    nameCStr = ccall(
        (:mongoc_collection_get_name, libmongoc),
        Ptr{UInt8}, (Ptr{Cvoid},),
        collection._wrap_
        )
    name = unsafe_string(nameCStr)
    print(io, "MongoCollection($name)")
end

baremodule MongoInsertFlags
    const None = 0
    const ContinueOnError = 1
end
export MongoInsertFlags

function insert(
    collection::MongoCollection,
    document::BSONObject;
    flags::Int = MongoInsertFlags.None
    )

    oid = get!(document, "_id", BSONOID())

    bsonError = BSONError()
    ccall(
        (:mongoc_collection_insert, libmongoc),
        Bool, (Ptr{Cvoid}, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{UInt8}),
        collection._wrap_,
        flags,
        document._wrap_,
        C_NULL,
        bsonError._wrap_
        ) || error("insert: $(string(bsonError))")

    return oid
end

insert(
    collection::MongoCollection,
    dict::Associative;
    flags::Int = MongoInsertFlags.None
    ) = insert(collection, BSONObject(dict), flags=flags)

insert(
    collection::MongoCollection,
    dict::NakedDict;
    flags::Int = MongoInsertFlags.None
    ) = insert(collection, BSONObject(Dict(dict)), flags=flags)
export insert

baremodule MongoUpdateFlags
    const None = 0
    const Upsert = 1
    const MultiUpdate = 2
end
export MongoUpdateFlags

function update(
    collection::MongoCollection,
    selector::BSONObject,
    change::BSONObject;
    flags::Int = MongoUpdateFlags.None
    )

    bsonError = BSONError()
    ccall(
        (:mongoc_collection_update, libmongoc),
        Bool, (Ptr{Cvoid}, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{UInt8}),
        collection._wrap_,
        flags,
        selector._wrap_,
        change._wrap_,
        C_NULL,
        bsonError._wrap_
        ) || error("update: $(string(bsonError))")
end

update(
    collection::MongoCollection,
    selector::Associative,
    change::Associative;
    flags::Int = MongoUpdateFlags.None
    ) = update(
        collection,
        BSONObject(selector),
        BSONObject(change),
        flags = flags
        )

update(c::MongoCollection, s::NakedDict, chg::NakedDict; kwargs...) =
    update(c, BSONObject(Dict(s)), BSONObject(Dict(chg)); kwargs...)
export update

baremodule MongoQueryFlags
    const None              = 0
    const TailableCursor    = 2
    const SlaveIO           = 4
    const OplogReplay       = 8
    const NoCursorTimeout   = 16
    const AwaitData         = 32
    const Exhaust           = 64
    const Partial           = 128
end
export MongoQueryFlags

function Base.find(
    collection::MongoCollection,
    selector::BSONObject,
    fields::BSONObject;
    flags::Int = MongoQueryFlags.None,
    skip::Int = 0,
    limit::Int = 0,
    batch_size::Int = 0
    )

    result = ccall(
        (:mongoc_collection_find, libmongoc),
        Ptr{Cvoid}, (Ptr{Cvoid}, Cint, UInt32, UInt32, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        collection._wrap_,
        flags,
        skip,
        limit,
        batch_size,
        selector._wrap_,
        fields._wrap_,
        C_NULL
        )
    result == C_NULL && error("mongoc_collection_find: failure")
    return MongoCursor( result )
end

Base.find(
    collection::MongoCollection,
    selector::Associative,
    fields::Associative;
    skip::Int = 0,
    limit::Int = 0,
    batch_size::Int = 0,
    flags::Int = MongoQueryFlags.None
    ) = find(
        collection,
        BSONObject(selector),
        BSONObject(fields),
        skip = skip,
        limit = limit,
        batch_size = batch_size,
        flags = flags
        )

function Base.find(
    collection::MongoCollection,
    selector::BSONObject;
    flags::Int = MongoQueryFlags.None,
    skip::Int = 0,
    limit::Int = 0,
    batch_size::Int = 0
    )

    result = ccall(
        (:mongoc_collection_find, libmongoc),
        Ptr{Cvoid}, (Ptr{Cvoid}, Cint, UInt32, UInt32, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        collection._wrap_,
        flags,
        skip,
        limit,
        batch_size,
        selector._wrap_,
        C_NULL,
        C_NULL
        )
    result == C_NULL && error("mongoc_collection_find: failure")
    return MongoCursor( result )
end

Base.find(
    collection::MongoCollection,
    selector::Associative;
    skip::Int = 0,
    limit::Int = 0,
    batch_size::Int = 0,
    flags::Int = MongoQueryFlags.None
    ) = find(
        collection,
        BSONObject(selector),
        skip = skip,
        limit = limit,
        batch_size = batch_size,
        flags = flags
        )

Base.find(c::MongoCollection, s::NakedDict; kwargs...) =
    find(c, BSONObject(Dict(s)); kwargs...)

Base.find(c::MongoCollection, s::NakedDict, f::NakedDict; kwargs...) =
    find(c, BSONObject(Dict(s)), BSONObject(Dict(f)); kwargs...)

function Base.count(
    collection::MongoCollection,
    queryBSON::BSONObject;
    skip::Int64 = 0,
    limit::Int64 = 0,
    flags::Int = MongoQueryFlags.None
    )

    bsonError = BSONError()
    result = ccall(
        (:mongoc_collection_count, libmongoc),
        Int64, (Ptr{Cvoid}, Cint, Ptr{Cvoid}, Int64, Int64, Ptr{Cvoid}, Ptr{UInt8}),
        collection._wrap_,
        flags,
        queryBSON._wrap_,
        skip,
        limit,
        C_NULL,
        bsonError._wrap_
        )
    result < 0 && error("count: $(string(bsonError))")
    return result
end

Base.count(
    collection::MongoCollection,
    query::Associative;
    skip::Int64 = 0,
    limit::Int64 = 0,
    flags::Int = MongoQueryFlags.None
    ) = count(
        collection,
        BSONObject(query),
        skip = skip,
        limit = limit,
        flags = flags
        )

Base.count(c::MongoCollection, s::NakedDict; kwargs...) =
    count(c, BSONObject(Dict(s)); kwargs...)

Base.count(c::MongoCollection) = count(c, BSONObject())

baremodule MongoDeleteFlags
    const None              = 0
    const SingleRemove      = 1
end
export MongoDeleteFlags

delete(
    collection::MongoCollection,
    selector::BSONObject;
    flags::Int = MongoDeleteFlags.None
    ) = begin
    bsonError = BSONError()
    result = ccall(
        (:mongoc_collection_delete, libmongoc),
        Bool, (Ptr{Cvoid}, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{UInt8}),
        collection._wrap_,
        flags,
        selector._wrap_,
        C_NULL,
        bsonError._wrap_
        )
    result < 0 && error("delete: $(string(bsonError))")
    return result
end

delete(
    collection::MongoCollection,
    selector::Associative;
    flags::Int = MongoDeleteFlags.None
    ) = delete(
        collection,
        BSONObject(selector),
        flags = flags
        )

delete(c::MongoCollection, s::NakedDict; kwargs...) =
    delete(c, BSONObject(Dict(s)); kwargs...)
export delete

destroy(collection::MongoCollection) =
    ccall(
        (:mongoc_collection_destroy, libmongoc),
        Cvoid, (Ptr{Cvoid},),
        collection._wrap_
        )


text:
using CloseOpenIntervals
using Documenter

DocMeta.setdocmeta!(CloseOpenIntervals, :DocTestSetup, :(using CloseOpenIntervals); recursive=true)

makedocs(;
    modules=[CloseOpenIntervals],
    authors="chriselrod <elrodc@gmail.com> and contributors",
    repo="https://github.com/JuliaSIMD/CloseOpenIntervals.jl/blob/{commit}{path}#{line}",
    sitename="CloseOpenIntervals.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaSIMD.github.io/CloseOpenIntervals.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaSIMD/CloseOpenIntervals.jl",
    devbranch = "main",
)


text:
module MakieGraphs

using AbstractPlotting
using GraphPlot
using LightGraphs
using Random

export NetworkPlot

@recipe(NetworkPlot, graph) do scene
    Attributes(
        nodecolor = :red,
        nodesize = 20,
        marker = Circle,
        edgecolor = :black,
        edgewidth = 2.5,
        seed = 39,
        c = 0.5
    )
end

function AbstractPlotting.plot!(np::NetworkPlot)
    graph = np[:graph]

    edges = Node(Point2f0[])
    nodes = Node(Point2f0[])

    function update_graph(new_graph, new_seed, new_c)
        empty!(edges[])
        empty!(nodes[])

        rng = MersenneTwister(new_seed)
        g = new_graph
        xpos, ypos = spring_layout(g, 2. * rand(rng, nv(g)) .- 1.0, 2 .* rand(rng, nv(g)) .- 1.0; C=new_c)

        for (x, y) in zip(xpos, ypos)
            push!(nodes[], Point2f0(x, -y))
        end
        for e in LightGraphs.edges(g)
            push!(edges[], nodes[][e.src])
            push!(edges[], nodes[][e.dst])
        end
        # trigger the changes
        edges[] = edges[]
        nodes[] = nodes[]
    end

    AbstractPlotting.Observables.onany(update_graph, graph, np.seed, np.c)

    update_graph(graph[], np.seed[], np.c[])

    linesegments!(np, edges,
                  color=np.edgecolor,
                  linewidth=np.edgewidth)
    scatter!(np, nodes,
             color=np.nodecolor,
             marker=np.marker,
             markersize=np.nodesize)
    np
end

end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Plots
using TaylorSeries

n = 15
x = range(-4π, 4π; length=1000)
plot(; legend=:outertopright)
plot!(x, sin.(x); label="sin(x)")
for k in 3:2:2n+1
    f(x) = evaluate(sin(Taylor1(k)), x)
    plot!(x, f.(x); label="k = $k", ls=:auto)
end
plot!(; ylim=(-10, 10), size=(600, 300))

# %%
using Plots
using QuadGK

function fourier_coeffs(f, N)
    a, b = eps(), 1-eps()
    c = Vector{Float64}(undef, 2N+1)
    c[1] = quadgk(f, a, b)[1]
    for k in 1:N
        c[k+1]   = quadgk(x -> cospi(2k*x)*f(x), a, b)[1]
        c[k+N+1] = quadgk(x -> sinpi(2k*x)*f(x), a, b)[1]
    end
    c
end

function fourier_sum(x, c, N)
    L = (length(c) - 1) ÷ 2
    c[1] + 2sum(c[k+1]*cospi(2k*x) + c[k+L+1]*sinpi(2k*x) for k in 1:N)
end

function fourier_plot(f, N=20)
    c = fourier_coeffs(f, N)
    x = range(-0.6, 1.6; length=1101)
    plot(; legend=:outertopright)
    f_mod1(x) = f(mod(x, 1))
    plot!(x, f_mod1.(x); label="f(mod(x, 1))")
    for k in (N, N÷5)
        f_k(x) = fourier_sum(x, c, k)
        plot!(x, f_k.(x); label="k = $k")
    end
    plot!(; size=(600, 300))
end

# %%
f(x) = x
fourier_plot(f)

# %%
f(x) = 0 ≤ x < 0.5 ? 1.0 : 0.0
fourier_plot(f)

# %%
f(x) = (x - 0.5)^2
fourier_plot(f)

# %%

max_stars_repo_path: 0005/Taylor and Fourier.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0005/Taylor and Fourier.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0005/Taylor and Fourier.jl

max_forks_repo_name: genkuroki/public


text:
###### update ####################

function update(ta::TimeArray{T,N,D}, tstamp::D, val::AbstractArray{T,N}) where {T,N,D}

    if length(ta) == 0
        throw(ArgumentError(
            "updating empty time arrays is not supported, " *
            "please use a scalable approach"))
    elseif tstamp < maximum(timestamp(ta))
        throw(ArgumentError("only appending operations supported"))
    else
        t    = vcat(timestamp(ta), tstamp)
        vals = vcat(values(ta), val)
        uta  = TimeArray(t, vals, colnames(ta), meta(ta))
    end
    uta
end

function update(ta::TimeArray{T,N,D}, tstamp::D, val::T) where {T,N,D}

    if length(ta) == 0
        throw(ArgumentError(
            "updating empty time arrays is not supported, " *
            "please use a scalable approach"))
    elseif tstamp < maximum(timestamp(ta))
        throw(ArgumentError("only appending operations supported"))
    else
        t    = vcat(timestamp(ta), tstamp)
        vals = vcat(values(ta), val)
        uta  = TimeArray(t, vals, colnames(ta), meta(ta))
    end
    uta
end

###### rename ####################

function rename(ta::TimeArray, colnames::Vector{Symbol})
    length(colnames) == size(ta, 2) || throw(DimensionMismatch("Colnames length mismatch"))
    TimeArray(ta, colnames = colnames, unchecked = true)
end

rename(ta::TimeArray, colnames::Symbol) = rename(ta, [colnames])


text:
module MetaGraphsNext

using JLD2
using Graphs

export MetaGraph, MetaDiGraph, MetaUndirectedGraph
export label_for, code_for, set_data
export weighttype, default_weight, weight_function
export MGFormat, DOTFormat

include("metagraph.jl")
include("metaundigraph.jl")
include("metadigraph.jl")
include("graphs.jl")
include("dict_utils.jl")
include("weights.jl")
include("persistence.jl")

end # module


text:
const FILE_FORMAT = ["2.2", "0", "8"]

immutable Physical
    num::Int
    dim::Int
    name::String
end

immutable Node
    x::Float64
    y::Float64
    z::Float64
end

immutable Element
    elmtype::Symbol
    phystag::Int
    geomtag::Int
    nodes::Vector{Int}
    function Element(elmtype::Int,phystag::Int,geomtag::Int,nodes::Vector{Int})
        if     elmtype==2; elmtype=:triangle
        elseif elmtype==3; elmtype=:quadrangle
        else error("The mesh can only contain 2D triangles or quadrangles!")
        end
        return new(elmtype,phystag,geomtag,nodes)
    end
end

immutable Mesh
    physicals::Vector{Physical}
    nodes::Vector{Node}
    elements::Vector{Element}
end

function get_shapes(m::Mesh)
    res = map(m.physicals) do p
        file = print_fastcap_file(m,p)
        return GmshShape(p.name,file)
    end
    return res
end

import Base.print
function print(io::IO,n::Node)
    return @printf(io,"% .6E % .6E % .6E",n.x,n.y,n.z)
end

function print_fastcap_file(m::Mesh,p::Physical)
    res = map(get_elms(m,p)) do elm
        nds = get_nodes(m,elm)
        if elm.elmtype==:triangle
            s = string("T 1 ",join(nds," ")) # join uses the print method above
        elseif elm.elmtype==:quadrangle
            s = string("Q 1 ",join(nds," "))
        else
            error("This shape is not supported by FastCap!")
        end
        return s
    end
    return join(["0 Gmsh physical: $(p.name)";res],NL)
end

function get_nodes(m::Mesh,e::Element)
    return m.nodes[e.nodes]
end

function get_elms(m::Mesh,p::Physical)
    return get_elms(m,p.num)
end

function get_elms(m::Mesh,p::Int)
    return collect(filter(e->e.phystag==p,m.elements))
end

function read_msh(filename::String)
    physicals,nodes,elms = open(parse_msh,filename,"r")
    return Mesh(physicals,nodes,elms)
end

function parse_msh(f::IOStream)
    local physicals,nodes,elms
    while !eof(f)
        header = strip( readline(f) )
        if header == "\$MeshFormat"
            fmt = split( readline(f) )
            @assert fmt == FILE_FORMAT
            @assert "\$EndMeshFormat" == strip(readline(f))
        elseif header == "\$PhysicalNames"
            nonames = parse(Int,readline(f))
            physicals = Vector{Physical}(nonames)
            for n = 1:nonames
                s = split(readline(f))
                @assert n == parse(Int,s[2])
                physicals[n] = Physical(n,
                    parse(Int,s[1]),
                    strip(s[3],'"'))
            end
            @assert "\$EndPhysicalNames" == strip(readline(f))
        elseif header == "\$Nodes"
            nonodes = parse(Int, readline(f))
            nodes = Vector{Node}(nonodes)
            for n = 1:nonodes
                s = split(readline(f))
                @assert n == parse(Int,s[1])
                nodes[n] = Node(parse(Float64,s[2]),parse(Float64,s[3]),parse(Float64,s[4]))
                #nodes[n] = Node(parse.([Float64],s[2:4])...)
            end
            @assert "\$EndNodes" == strip(readline(f))
        elseif header == "\$Elements"
            noelms = parse(Int, readline(f))
            elms = Vector{Element}(noelms)
            for n = 1:noelms
                s = [parse(Int,x) for x in split(readline(f))]
                #s = parse.([Int],split(readline(f)))
                @assert n == s[1]
                notags = s[3]
                elms[n] = Element(s[2],s[4],s[5],s[6:end])
            end
            @assert "\$EndElements" == strip(readline(f))
        elseif header == "\$Comment"
            while true
                line = readline(f)
                if strip(line) == "\$EndComment"
                    break
                end
            end
        else
            msg = string("Unknown section header: ", line)
            error(msg)
        end
    end
    return physicals,nodes,elms
end

max_stars_repo_name: jebej/FastCap

max_issues_repo_name: jebej/FastCap

max_forks_repo_name: jebej/FastCap


text:
module ImageMetadata

# The order here is designed to avoid an ambiguity warning in convert,
# see the top of ImageAxes
using ImageAxes
using ImageCore, Colors, FixedPointNumbers
using ColorVectorSpace   # for overriding math operations with Gray/RGB
import AxisArrays

import Base: +, -, *, /
import Base: permutedims

const ViewIndex = Union{Base.ViewIndex, Colon}

export
    # types
    ImageMeta,

    # functions
    copyproperties,
    data,
    properties,
    shareproperties,
    spatialproperties

#### types and constructors ####

# Concrete types
"""
`ImageMeta` is an AbstractArray that can have metadata, stored in a dictionary.

Construct an image with `ImageMeta(A, props)` (for a properties dictionary
`props`), or with `ImageMeta(A, prop1=val1, prop2=val2, ...)`.
"""
mutable struct ImageMeta{T,N,A<:AbstractArray,P<:AbstractDict{String,Any}} <: AbstractArray{T,N}
    data::A
    properties::P

    function ImageMeta{T,N,A,P}(data::AbstractArray, properties::P) where {T,N,A,P}
        new{T,N,A,P}(data, properties)
    end
end
ImageMeta(data::AbstractArray{T,N}, props::AbstractDict{String,Any}) where {T,N} = ImageMeta{T,N,typeof(data),typeof(props)}(data,props)
ImageMeta(data::AbstractArray{T,N}, props::Dict{String,Any}) where {T,N} = ImageMeta{T,N,typeof(data),typeof(props)}(data,props)
ImageMeta(data::AbstractArray, props::Dict{<:AbstractString}) = ImageMeta(data, convert(Dict{String,Any}, props))
ImageMeta(data::AbstractArray; kwargs...) = ImageMeta(data, kwargs2dict(kwargs))
ImageMeta(data::AbstractArray, props::AbstractDict) = throw(ArgumentError("properties must be an AbstractDict{String,Any}"))

const ImageMetaArray{T,N,A<:Array} = ImageMeta{T,N,A}
const ImageMetaAxis{T,N,A<:AxisArray} = ImageMeta{T,N,A}

Base.size(A::ImageMeta) = size(A.data)

datatype(::Type{ImageMeta{T,N,A,P}}) where {T,N,A<:AbstractArray,P} = A

Base.IndexStyle(::Type{M}) where {M<:ImageMeta} = IndexStyle(datatype(M))

AxisArrays.HasAxes(A::ImageMetaAxis) = AxisArrays.HasAxes{true}()

# getindex and setindex!
for AType in (ImageMeta, ImageMetaAxis)
    @eval begin
        @inline function Base.getindex(img::$AType{T,1}, i::Int) where T
            @boundscheck checkbounds(img.data, i)
            @inbounds ret = img.data[i]
            ret
        end
        @inline function Base.getindex(img::$AType, i::Int)
            @boundscheck checkbounds(img.data, i)
            @inbounds ret = img.data[i]
            ret
        end
        @inline function Base.getindex(img::$AType{T,N}, I::Vararg{Int,N}) where {T,N}
            @boundscheck checkbounds(img.data, I...)
            @inbounds ret = img.data[I...]
            ret
        end

        @inline function Base.setindex!(img::$AType{T,1}, val, i::Int) where T
            @boundscheck checkbounds(img.data, i)
            @inbounds img.data[i] = val
            val
        end
        @inline function Base.setindex!(img::$AType, val, i::Int)
            @boundscheck checkbounds(img.data, i)
            @inbounds img.data[i] = val
            val
        end
        @inline function Base.setindex!(img::$AType{T,N}, val, I::Vararg{Int,N}) where {T,N}
            @boundscheck checkbounds(img.data, I...)
            @inbounds img.data[I...] = val
            val
        end
    end
end

@inline function Base.getindex(img::ImageMetaAxis, ax::Axis, I...)
    result = img.data[ax, I...]
    maybe_wrap(img, result)
end
@inline function Base.getindex(img::ImageMetaAxis, i::Union{Integer,AbstractVector,Colon}, I...)
    result = img.data[i, I...]
    maybe_wrap(img, result)
end
maybe_wrap(img::ImageMeta{T}, result::T) where T = result
maybe_wrap(img::ImageMeta{T}, result::AbstractArray{T}) where T = copyproperties(img, result)

@inline function Base.setindex!(img::ImageMetaAxis, val, ax::Axis, I...)
    setindex!(img.data, val, ax, I...)
end
@inline function Base.setindex!(img::ImageMetaAxis, val, i::Union{Integer,AbstractVector,Colon}, I...)
    setindex!(img.data, val, i, I...)
end

Base.view(img::ImageMeta, ax::Axis, I...) = shareproperties(img, view(img.data, ax, I...))
Base.view(img::ImageMeta{T,N}, I::Vararg{ViewIndex,N}) where {T,N} = shareproperties(img, view(img.data, I...))
Base.view(img::ImageMeta, i::ViewIndex) = shareproperties(img, view(img.data, i))
Base.view(img::ImageMeta, I::Vararg{ViewIndex,N}) where {N} = shareproperties(img, view(img.data, I...))

Base.getindex(img::ImageMeta, propname::AbstractString) = img.properties[propname]

Base.setindex!(img::ImageMeta, X, propname::AbstractString) = setindex!(img.properties, X, propname)

Base.copy(img::ImageMeta) = ImageMeta(copy(img.data), deepcopy(img.properties))

Base.convert(::Type{ImageMeta}, A::ImageMeta) = A
Base.convert(::Type{ImageMeta}, A::AbstractArray) = ImageMeta(A)
Base.convert(::Type{ImageMeta{T}}, A::ImageMeta{T}) where {T} = A
Base.convert(::Type{ImageMeta{T}}, A::ImageMeta) where {T} = shareproperties(A, convert(Array{T}, A.data))
Base.convert(::Type{ImageMeta{T}}, A::AbstractArray) where {T} = ImageMeta(convert(Array{T}, A))

# copy properties
function Base.copy!(imgdest::ImageMeta, imgsrc::ImageMeta, prop1::AbstractString, props::AbstractString...)
    imgdest[prop1] = deepcopy(imgsrc[prop1])
    for p in props
        imgdest[p] = deepcopy(imgsrc[p])
    end
    imgdest
end

# similar
Base.similar(img::ImageMeta, ::Type{T}, shape::Dims) where {T} = ImageMeta(similar(img.data, T, shape), deepcopy(img.properties))
Base.similar(img::ImageMetaAxis, ::Type{T}) where {T} = ImageMeta(similar(img.data, T), deepcopy(img.properties))

"""
    copyproperties(img::ImageMeta, data) -> imgnew

Create a new "image," copying the properties dictionary of `img` but
using the data of the AbstractArray `data`. Note that changing the
properties of `imgnew` does not affect the properties of `img`.

See also: [`shareproperties`](@ref).
"""
copyproperties(img::ImageMeta, data::AbstractArray) =
    ImageMeta(data, deepcopy(img.properties))

"""
    shareproperties(img::ImageMeta, data) -> imgnew

Create a new "image," reusing the properties dictionary of `img` but
using the data of the AbstractArray `data`. The two images have
synchronized properties; modifying one also affects the other.

See also: [`copyproperties`](@ref).
"""
shareproperties(img::ImageMeta, data::AbstractArray) = ImageMeta(data, img.properties)

# Delete a property!
Base.delete!(img::ImageMeta, propname::AbstractString) = delete!(img.properties, propname)

# Iteration
# Defer to the array object in case it has special iteration defined
Base.iterate(img::ImageMeta) = Base.iterate(data(img))
Base.iterate(img::ImageMeta, s) = Base.iterate(data(img), s)

# Show
const emptyset = Set()
function showim(io::IO, img::ImageMeta)
    IT = typeof(img)
    print(io, eltype(img).name.name, " ImageMeta with:\n  data: ", summary(img.data), "\n  properties:")
    showdictlines(io, img.properties, get(img, "suppress", emptyset))
end
Base.show(io::IO, img::ImageMeta) = showim(io, img)
Base.show(io::IO, ::MIME"text/plain", img::ImageMeta) = showim(io, img)

function Base.reinterpret(::Type{T}, img::ImageMeta) where {T}
    shareproperties(img, reinterpret(T, img.data))
end

"""
    data(img::ImageMeta) -> array

Extract the data from `img`, omitting the properties
dictionary. `array` shares storage with `img`, so changes to one
affect the other.

See also: [`properties`](@ref).
"""
ImageAxes.data(img::ImageMeta) = img.data   # fixme when deprecation is removed from ImageCore

function ImageCore.permuteddimsview(A::ImageMeta, perm)
    ip = sortperm([perm...][[coords_spatial(A)...]])  # the inverse spatial permutation
    permutedims_props!(copyproperties(A, permuteddimsview(A.data, perm)), ip)
end
ImageCore.channelview(A::ImageMeta) = shareproperties(A, channelview(A.data))
ImageCore.colorview(::Type{C}, A::ImageMeta{T,N}) where {C<:Colorant,T,N} = shareproperties(A, colorview(C, A.data))
ImageCore.colorview(::Type{ARGB32}, A::ImageMeta{T,N}) where {T,N} = shareproperties(A, colorview(ARGB32, A.data))
ImageCore.rawview(A::ImageMeta{T}) where {T<:Real} = shareproperties(A, rawview(A.data))
ImageCore.normedview(::Type{T}, A::ImageMeta{S}) where {T<:FixedPoint,S<:Unsigned} = shareproperties(A, normedview(T, A.data))

# AxisArrays functions
AxisArrays.axes(img::ImageMetaAxis) = AxisArrays.axes(img.data)
AxisArrays.axisdim(img::ImageMetaAxis, ax) = axisdim(img.data, ax)
AxisArrays.axisnames(img::ImageMetaAxis) = axisnames(img.data)
AxisArrays.axisvalues(img::ImageMetaAxis) = axisvalues(img.data)

#### Properties ####

"""
    properties(imgmeta) -> props

Extract the properties dictionary `props` for `imgmeta`. `props`
shares storage with `img`, so changes to one affect the other.

See also: [`data`](@ref).
"""
properties(img::ImageMeta) = img.properties

Base.haskey(img::ImageMeta, k::AbstractString) = haskey(img.properties, k)

Base.get(img::ImageMeta, k::AbstractString, default) = get(img.properties, k, default)

# So that defaults don't have to be evaluated unless they are needed,
# we also define a @get macro (thanks Toivo Hennington):
struct IMNothing end   # to avoid confusion in the case where dict[key] === nothing
macro get(img, k, default)
    quote
        img, k = $(esc(img)), $(esc(k))
        val = get(img.properties, k, IMNothing())
        return isa(val, IMNothing) ? $(esc(default)) : val
    end
end

ImageAxes.timeaxis(img::ImageMetaAxis) = timeaxis(data(img))
ImageAxes.timedim(img::ImageMetaAxis) = timedim(data(img))
ImageAxes.colordim(img::ImageMetaAxis) = colordim(data(img))

ImageCore.pixelspacing(img::ImageMeta) = pixelspacing(data(img))

"""
    spacedirections(img)

Using ImageMetadata, you can set this property manually. For example, you
could indicate that a photograph was taken with the camera tilted
30-degree relative to vertical using

```
img["spacedirections"] = ((0.866025,-0.5),(0.5,0.866025))
```

If not specified, it will be computed from `pixelspacing(img)`, placing the
spacing along the "diagonal".  If desired, you can set this property in terms of
physical units, and each axis can have distinct units.
"""
ImageCore.spacedirections(img::ImageMeta) = @get img "spacedirections" spacedirections(data(img))

ImageCore.sdims(img::ImageMetaAxis) = sdims(data(img))

ImageCore.coords_spatial(img::ImageMetaAxis) = coords_spatial(data(img))

ImageCore.spatialorder(img::ImageMetaAxis) = spatialorder(data(img))

ImageAxes.nimages(img::ImageMetaAxis) = nimages(data(img))

ImageCore.size_spatial(img::ImageMetaAxis) = size_spatial(data(img))

ImageCore.indices_spatial(img::ImageMetaAxis) = indices_spatial(data(img))

ImageCore.assert_timedim_last(img::ImageMetaAxis) = assert_timedim_last(data(img))

#### Permutations over dimensions ####

"""
    permutedims(img, perm, [spatialprops])

When permuting the dimensions of an ImageMeta, you can optionally
specify that certain properties are spatial and they will also be
permuted. `spatialprops` defaults to `spatialproperties(img)`.
"""
permutedims

function permutedims_props!(ret::ImageMeta, ip, spatialprops=spatialproperties(ret))
    if !isempty(spatialprops)
        for prop in spatialprops
            if haskey(ret, prop)
                a = ret.properties[prop]
                if isa(a, AbstractVector)
                    ret.properties[prop] = a[ip]
                elseif isa(a, Tuple)
                    ret.properties[prop] = a[ip]
                elseif isa(a, AbstractMatrix) && size(a,1) == size(a,2)
                    ret.properties[prop] = a[ip,ip]
                else
                    error("Do not know how to handle property ", prop)
                end
            end
        end
    end
    ret
end

function permutedims(img::ImageMetaAxis, perm)
    p = AxisArrays.permutation(perm, axisnames(img.data))
    ip = sortperm([p...][[coords_spatial(img)...]])
    permutedims_props!(copyproperties(img, permutedims(img.data, p)), ip)
end
function permutedims(img::ImageMeta, perm)
    ip = sortperm([perm...][[coords_spatial(img)...]])
    permutedims_props!(copyproperties(img, permutedims(img.data, perm)), ip)
end

# Note: `adjoint` does not recurse into ImageMeta properties.
function Base.adjoint(img::ImageMeta{T,2}) where {T<:Real}
    ip = sortperm([2,1][[coords_spatial(img)...]])
    permutedims_props!(copyproperties(img, adjoint(img.data)), ip)
end

function Base.adjoint(img::ImageMeta{T,1}) where T<:Real
    check_empty_spatialproperties(img)
    copyproperties(img, img.data')
end

"""
    spatialproperties(img)

Return a vector of strings, containing the names of properties that
have been declared "spatial" and hence should be permuted when calling
`permutedims`.  Declare such properties like this:

    img["spatialproperties"] = ["spacedirections"]
"""
spatialproperties(img::ImageMeta) = @get img "spatialproperties" ["spacedirections"]

function check_empty_spatialproperties(img)
    sp = spatialproperties(img)
    for prop in sp
        if haskey(img, prop)
            error("spatialproperties must be empty, have $prop")
        end
    end
    nothing
end

#### Low-level utilities ####
function showdictlines(io::IO, dict::AbstractDict, suppress::Set)
    for (k, v) in dict
        if k == "suppress"
            continue
        end
        if !in(k, suppress)
            print(io, "\n    ", k, ": ")
            print(IOContext(io, :compact => true), v)
        else
            print(io, "\n    ", k, ": <suppressed>")
        end
    end
end
showdictlines(io::IO, dict::AbstractDict, prop::String) = showdictlines(io, dict, Set([prop]))

# printdictval(io::IO, v) = print(io, v)
# function printdictval(io::IO, v::Vector)
#     for i = 1:length(v)
#         print(io, " ", v[i])
#     end
# end

# converts keyword argument to a dictionary
function kwargs2dict(kwargs)
    d = Dict{String,Any}()
    for (k, v) in kwargs
        d[string(k)] = v
    end
    return d
end

include("operators.jl")

end # module


text:
# ---
# title: 599. Minimum Index Sum of Two Lists
# id: problem599
# author: zhwang
# date: 2022-03-14
# difficulty: Easy
# categories: Hash Table
# link: <https://leetcode.com/problems/minimum-index-sum-of-two-lists/description/>
# hidden: true
# ---
#
# Suppose Andy and Doris want to choose a restaurant for dinner, and they both
# have a list of favorite restaurants represented by strings.
#
# You need to help them find out their **common interest** with the **least list
# index sum**. If there is a choice tie between answers, output all of them with
# no order requirement. You could assume there always exists an answer.
#
#
#
# **Example 1:**
#
#
#
#     Input: list1 = ["Shogun","Tapioca Express","Burger King","KFC"], list2 = ["Piatti","The Grill at Torrey Pines","Hungry Hunter Steakhouse","Shogun"]
#     Output: ["Shogun"]
#     Explanation: The only restaurant they both like is "Shogun".
#
#
# **Example 2:**
#
#
#
#     Input: list1 = ["Shogun","Tapioca Express","Burger King","KFC"], list2 = ["KFC","Shogun","Burger King"]
#     Output: ["Shogun"]
#     Explanation: The restaurant they both like and have the least index sum is "Shogun" with index sum 1 (0+1).
#
#
# **Example 3:**
#
#
#
#     Input: list1 = ["Shogun","Tapioca Express","Burger King","KFC"], list2 = ["KFC","Burger King","Tapioca Express","Shogun"]
#     Output: ["KFC","Burger King","Tapioca Express","Shogun"]
#
#
# **Example 4:**
#
#
#
#     Input: list1 = ["Shogun","Tapioca Express","Burger King","KFC"], list2 = ["KNN","KFC","Burger King","Tapioca Express","Shogun"]
#     Output: ["KFC","Burger King","Tapioca Express","Shogun"]
#
#
# **Example 5:**
#
#
#
#     Input: list1 = ["KFC"], list2 = ["KFC"]
#     Output: ["KFC"]
#
#
#
#
# **Constraints:**
#
#   * `1 <= list1.length, list2.length <= 1000`
#   * `1 <= list1[i].length, list2[i].length <= 30`
#   * `list1[i]` and `list2[i]` consist of spaces `' '` and English letters.
#   * All the stings of `list1` are **unique**.
#   * All the stings of `list2` are **unique**.
#
#
## @lc code=start
using LeetCode

function find_restaurant(list1::Vector{String}, list2::Vector{String})::Vector{String}
    dic = Dict(st => (i - 1) for (i, st) in enumerate(list1))
    cost, res = typemax(1), String[]
    for (i, s) in enumerate(list2)
        !haskey(dic, s) && continue
        newcost = i - 1 + dic[s]
        if newcost < cost
            cost, res = newcost, [s]
        elseif newcost == cost
            push!(res, s)
        end
    end
    return res
end

## @lc code=end


text:
function TT_SVD(tensor::Array{Float64},ϵ::Float64)
########################################################################
#   Computes the cores of a TT for the given tensor and accuracy (acc)
#   Resources:
#   V. Oseledets: Tensor-Train Decomposition, 2011, p.2301: Algorithm 1
#   April 2021, Clara Menzen
########################################################################
    N = ndims(tensor);
    cores = Vector{Array{Float64,3}}(undef,N);
    frobnorm = norm(tensor);

    δ = ϵ / sqrt(N-1) * frobnorm;
    err2 = 0;
    rprev = 1;
    sizes = size(tensor);
    C = reshape( tensor, (sizes[1], Int(length(tensor) / sizes[1]) ));
    for k = 1 : N-1
        # truncated svd
        F   = svd!(C);
        rcurr = length(F.S);

        sv2 = cumsum(reverse(F.S).^2);
        tr  = Int(findfirst(sv2 .> δ^2))-1;
        if tr > 0
            rcurr = length(F.S) - tr;
            err2 += sv2[tr];
        end

        # new core
        cores[k] = reshape(F.U[:,1:rcurr],(rprev,sizes[k],rcurr));
        rprev    = rcurr;
        C        = Diagonal(F.S[1:rcurr])*F.Vt[1:rcurr,:];
        C        = reshape(C,(rcurr*sizes[k+1], Int(length(C) / (rcurr*sizes[k+1])) ) );
    end
    cores[N] = reshape(C,(rprev,sizes[N],1));
    return MPT(cores), sqrt(err2)/frobnorm
end

function TT_SVD(tensor::Array{Float64},d::Int,ϵ::Float64)
    ########################################################################
    # for a symmetric matrix
        cores    = Vector{Array{Float64,3}}(undef,d);
        S        = zeros(2)
        frobnorm = norm(tensor);

        δ     = ϵ / sqrt(d-1) * frobnorm;
        err2  = 0;
        rprev = 1;
        sizes = size(tensor);
        C     = reshape( tensor, (sizes[1], Int(length(tensor) / sizes[1]) ));
        for k = 1 : d
            # truncated svd
            F   = svd!(C);
            rcurr = length(F.S);

            sv2 = cumsum(reverse(F.S).^2);
            tr  = Int(findfirst(sv2 .> δ^2))-1;
            if tr > 0
                rcurr = length(F.S) - tr;
                err2 += sv2[tr];
            end

            # new core

            cores[k] = reshape(F.U[:,1:rcurr],(rprev,sizes[k],rcurr));
            rprev    = rcurr;
            C        = Diagonal(F.S[1:rcurr])*F.Vt[1:rcurr,:];
            C        = reshape(C,(rcurr*sizes[k+1], Int(length(C) / (rcurr*sizes[k+1])) ) );
            S        = F.S[1:rcurr]
        end
        return MPT(cores), S, sqrt(err2)/frobnorm
    end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
using TSne  #to install run command Pkg.clone("git://github.com/lejon/TSne.jl.git")
using DataFrames

# This is from https://github.com/lejon/TSne.jl/blob/35bad1c5cfd78d44e5d89bbca89a85555943b9db/examples/demo-csv.jl

# USAGE: julia bin/create_tsne_plot.jl haveheader --labelcolname=pathway_label /home/awright/git/PGM/gecco/mpbiopath_results.summary.for_tsne.tsv

doc = """Use t-SNE to generate a PDF called myplot.pdf from an input CSV file. Default assumption is to have no header and no labels. If these are available in the CSV these must be given as arguments.
Usage:
  demo-csv.jl <filename>
  demo-csv.jl [--labelcol=<col>] <filename>
  demo-csv.jl haveheader [--labelcol=<col>] <filename>
  demo-csv.jl [--labelcolname=<colname>] <filename>
  demo-csv.jl haveheader [--labelcolname=<colname>] <filename>
Options:
  -h --help     Show this screen.
  --version     Show version.
  --filename=   Path to CSV file
  --noheader    The CSV file does not have a header row
  --nolabel    The CSV file does not have a label column
"""

function normalize(A)
	for col in 1:size(A)[2]
        	std(A[:,col]) == 0 && continue
        	A[:,col] = (A[:,col]-mean(A[:,col])) / std(A[:,col])
	end
	A
end

using DocOpt

arguments = docopt(doc, version=v"0.0.1")
dump(arguments)

if nothing==arguments["--labelcol"]
	lblcol = -1
else
	lblcol = parse(Int64,arguments["--labelcol"])
end

df = readtable(arguments["<filename>"], header = nothing!=arguments["haveheader"], separator = '\t')

if nothing!=arguments["--labelcolname"]
	lblcol = find(x -> x==symbol(arguments["--labelcolname"]),names(df))[1]
end

println("Data is $df")
if lblcol>0
	labels = df[:,lblcol]
end

dataset = df[filter(x -> x!=lblcol,1:ncol(df)),]

data = float(convert(Array,dataset))
# Normalize the data, this should be done if there are large scale differences in the dataset
X = normalize(data)

# Run t-SNE
Y = tsne(X, 2, 50, 1000, 20.0)

using Gadfly
if lblcol>0
	theplot = plot(x=Y[:,1], y=Y[:,2], color=labels)
else
	theplot = plot(x=Y[:,1], y=Y[:,2])
end
draw(PDF("tsne-plot2.pdf", 8inch, 6inch), theplot)

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
include("../src/DataPackage.jl")
using DataPackage
import DataPackage: Package, read

data_url = "https://datahub.io/core/pharmaceutical-drug-spending/datapackage.json"

# to load Data Package into storage
package = Package(data_url)

# to load only tabular data
resources = package.resources
for resource in resources
    if resource.profile == "tabular-data-resource"
        data = read(resource)
        println(data)
    end
end


text:
"""Sieve of Eratosthenes function docstring"""
function es(n::Int) # one integer argument
    isprime = trues(n) # n-elt vector of true
    isprime[1] = false # since 1 is not prime
    for i in 2:isqrt(n) # loop odd ints less or equal to sqrt(n)
        if isprime[i] # conditonal evaluation
            for j in i^2:i:n # sequence with step i
                isprime[j] = false
            end
        end
    end
    return filter(x -> isprime[x], 1:n) # filter using anonymous fn
end

println(es(100)) # print all primes less or equal to 100
@time length(es(10^6)) # check fn execution time and memory usage

max_stars_repo_name: mxxo/numerics

max_issues_repo_name: mxxo/numerics

max_forks_repo_name: mxxo/numerics


text:
module ProjectiveGeometricAlgebra2d
#https://bivector.net/2DPGA.pdf
using StaticArrays
using ArgCheck
using Random

struct PGA2d{T}
    e::T
    e0::T
    e1::T
    e2::T
    e12::T
    e02::T
    e01::T
    e012::T
end

function Base.show(io::IO, o::PGA2d)
    print(io, "PGA2d")
    print(io, "(")
    pieces = String[]
    for pname in propertynames(o)
        val = getproperty(o, pname)
        if val != 0
            push!(pieces, "$pname=$val")
        end
    end
    print(io, join(pieces, ", "))
    print(io, ")")
end

function _serialize(x::PGA2d)
    @SVector[
        x.e,x.e0,x.e1,x.e2,x.e12,x.e02,x.e01,x.e012,
    ]
end

function Base.isapprox(x::PGA2d, y::PGA2d; kw...)
    isapprox(_serialize(x), _serialize(y); kw...)
end

@inline function PGA2d(;e=0, e0=0, e1=0, e2=0, e12=0, e02=0, e01=0, e012=0)
    e, e0, e1, e2, e12, e02, e01, e012 = promote(e, e0, e1, e2, e12, e02, e01, e012)
    PGA2d(
        e   ,
        e0  ,
        e1  ,
        e2  ,
        e12 ,
        e02 ,
        e01 ,
        e012,
    )
end

function wedge(x,y)
    PGA2d(
        e = x.e*y.e,

        e0 = x.e*y.e0 + x.e0*y.e,
        e1 = x.e*y.e1 + x.e1*y.e,
        e2 = x.e*y.e2 + x.e2*y.e,

        e12 = x.e*y.e12 + x.e1*y.e2 - x.e2*y.e1 + x.e12*y.e,
        e02 = x.e*y.e02 + x.e0*y.e2 - x.e2*y.e0 + x.e02*y.e,
        e01 = x.e*y.e01 + x.e0*y.e1 - x.e1*y.e0 + x.e01*y.e,

        e012 = x.e*y.e012 +
        x.e0*y.e12 - x.e1*y.e02 + x.e2*y.e01 +
        x.e12*y.e0 - x.e02*y.e1 + x.e01*y.e2 +
        x.e012 * y.e,
    )
end

function geomul(x,y)
    PGA2d(
        e = x.e*y.e + x.e1*y.e1 + x.e2*y.e2 - x.e12*y.e12,

        e0 = x.e*y.e0 + x.e0*y.e - x.e1*y.e01 - x.e2*y.e02 - x.e12*y.e012 + x.e02*y.e2 + x.e01*y.e1 - x.e012*y.e12,
        e1 = x.e*y.e1 + x.e1*y.e - x.e2*y.e12 + x.e12*y.e2,
        e2 = x.e*y.e2 + x.e2*y.e + x.e1*y.e12 - x.e12*y.e1,

        e12 = x.e*y.e12 + x.e1*y.e2 - x.e2*y.e1 + x.e12*y.e,
        e02 = x.e*y.e02 + x.e0*y.e2 - x.e1*y.e012 - x.e2*y.e0 - x.e12*y.e01 + x.e02*y.e + x.e01*y.e12 - x.e012*y.e1,
        e01 = x.e*y.e01 + x.e0*y.e1 - x.e1*y.e0 + x.e2*y.e012 + x.e12*y.e02 - x.e02*y.e12 + x.e01*y.e + x.e012*y.e2,
        e012 = x.e*y.e012 + x.e0*y.e12 - x.e1*y.e02 + x.e2*y.e01 + x.e12*y.e0 - x.e02*y.e1 + x.e01*y.e2 + x.e012*y.e,
    )
end

@inline function dual(x::PGA2d)
    PGA2d(
        e   =  x.e012,
        e0  =  x.e12,
        e1  = -x.e02,
        e2  =  x.e01,
        e12 =  x.e0,
        e02 = -x.e1,
        e01 =  x.e2,
        e012= x.e,
    )
end

function reverse(x::PGA2d)
    PGA2d(
        e    = e    ,

        e0   = e0   ,
        e1   = e1   ,
        e2   = e2   ,

        e12  = -e12 ,
        e02  = -e02 ,
        e01  = -e01 ,

        e012 = -e012,
    )

end

function pga_point(v)
    @argcheck length(v) == 2
    x,y = v
    PGA2d(e02=-x, e01=y, e12=1)
end

function pga_line(line)
    # ax + by + c = 0
    PGA2d(e1=line.a, e2=line.b, e0=line.c)
end

function euclidean_line_abc(pga)
    (a=pga.e1, b=pga.e2, c=pga.e0)
end

function euclidean_point(pga)
    w = pga.e12
    x = - pga.e02 / w
    y = pga.e01 / w
    @SVector[x,y]
end

function join_points(pt1, pt2)
    l_pga = anti_wedge(pga_point(pt1), pga_point(pt2))
    euclidean_line_abc(l_pga)
end

function meet_lines(l1, l2)
    pt_pga = wedge(pga_line(l1), pga_line(l2))
    return euclidean_point(pt_pga)
end

function anti_wedge_using_dual(x1,x2)
    dual(wedge(dual(x1), dual(x2)))
end

function randn_pga(rng=Random.GLOBAL_RNG)
    blades = ntuple(_->randn(rng), Val(8))
    PGA2d(blades...)
end

function rand_pga(args...)
    blades = ntuple(_->rand(args...), Val(8))
    PGA2d(blades...)
end

function anti_wedge(x1,x2)
    # TODO more efficient implementation
    anti_wedge_using_dual(x1,x2)
end

function pseudo_one(T=Int)
    PGA2d(e012=one(T))
end

end#module