text:
import Base: @__doc__
using PowerDynamics: @DynamicNode
import PowerDynamics:
    construct_vertex, dimension, symbolsof, showdefinition, AbstractNode
using LinearAlgebra: Diagonal
using NetworkDynamics: ODEVertex

begin
    @__doc__ struct DGUnitPLL <: AbstractNode
        I_r
        K_pll
        Y_n
    end
    DGUnitPLL(; I_r, K_pll, Y_n) = DGUnitPLL(I_r, K_pll, Y_n)
    function construct_vertex(par::DGUnitPLL)
        I_r = par.I_r
        K_pll = par.K_pll
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            # current from incoming/outgoing edges, nodal shunt
            i = total_current(e_s, e_d) + Y_n * u
            θ = x[3]

            # PLL
            # (x[2] * cos(θ) - x[1] * sin(θ))
            dθ = -K_pll * abs(u) * sin(θ - angle(u))

            du = i - I_r * exp(θ * 1im) # = 0 constraint

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                dx[3] = dθ
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(
            f! = rhs!,
            dim = 3,
            mass_matrix = Diagonal([0, 0, 1]),
            sym = Symbol[:u_r, :u_i, :θ],
        )
    end
    symbolsof(::DGUnitPLL) = begin
        [:u_r, :u_i, :θ]
    end
    dimension(::DGUnitPLL) = begin
        3
    end
end

begin
    @__doc__ struct DGUnitPLLPQ <: AbstractNode
        I_r
        K_pll
        S_pq
        Y_n
    end
    DGUnitPLLPQ(; I_r, K_pll, S_pq, Y_n) = DGUnitPLLPQ(I_r, K_pll, S_pq, Y_n)
    function construct_vertex(par::DGUnitPLLPQ)
        I_r = par.I_r
        K_pll = par.K_pll
        S_pq = par.S_pq
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            # current from incoming/outgoing edges, nodal shunt, PQ load background
            i = total_current(e_s, e_d) + Y_n * u - conj(S_pq) / conj(u)
            θ = x[3]

            # PLL
            # (x[2] * cos(θ) - x[1] * sin(θ))
            dθ = -K_pll * abs(u) * sin(θ - angle(u))

            du = i - I_r * exp(θ * 1im) # = 0 constraint

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                dx[3] = dθ
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(
            f! = rhs!,
            dim = 3,
            mass_matrix = Diagonal([0, 0, 1]),
            sym = Symbol[:u_r, :u_i, :θ],
        )
    end
    symbolsof(::DGUnitPLLPQ) = begin
        [:u_r, :u_i, :θ]
    end
    dimension(::DGUnitPLLPQ) = begin
        3
    end
end

begin
    @__doc__ struct DGUnitPLLPQTracking <: AbstractNode
        K_pll
        K_PT1
        T_PT1
        S_pq
        Y_n
    end
    DGUnitPLLPQTracking(; K_pll, K_PT1, T_PT1, S_pq, Y_n) =
        DGUnitPLLPQTracking(K_pll, K_PT1, T_PT1, S_pq, Y_n)
    function construct_vertex(par::DGUnitPLLPQTracking)
        K_pll = par.K_pll
        K_PT1 = par.K_PT1
        T_PT1 = par.T_PT1
        S_pq = par.S_pq
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            # current from incoming/outgoing edges, nodal shunt, PQ load background
            i = total_current(e_s, e_d) + Y_n * u - conj(S_pq) / conj(u)
            Vamp = abs(u)

            θ = x[3]
            P_int = x[4]
            Q_int = x[5]
            P_g = x[6]
            Q_g = x[7]

            # power flow tracking
            ΔP, ΔQ = p

            dP_int = -ΔP
            dQ_int = -ΔQ
            dP_g = (K_PT1 * P_int / 2.0 - P_g) / T_PT1
            dQ_g = (K_PT1 * Q_int / 2.0 - Q_g) / T_PT1

            # current set point adjustment
            I_r = complex(P_g, -Q_g) / Vamp

            # PLL
            # (x[2] * cos(θ) - x[1] * sin(θ))
            dθ = -K_pll * Vamp * sin(θ - angle(u))

            du = i - I_r * exp(θ * 1im) # = 0 constraint

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                dx[3] = dθ
                dx[4] = dP_int
                dx[5] = dQ_int
                dx[6] = dP_g
                dx[7] = dQ_g
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(
            f! = rhs!,
            dim = 7,
            mass_matrix = Diagonal([0, 0, 1, 1, 1, 1, 1]),
            sym = Symbol[:u_r, :u_i, :θ, :P_int, :Q_int, :P_g, :Q_g],
        )
    end
    symbolsof(::DGUnitPLLPQTracking) = begin
        [:u_r, :u_i, :θ, :P_int, :Q_int, :P_g, :Q_g]
    end
    dimension(::DGUnitPLLPQTracking) = begin
        7
    end
end

begin
    @__doc__ struct DGUnitTrackingFRT <: AbstractNode
        K_pll
        K_PT1
        T_PT1
        K_FRT
        I_max
        Vref
        Vdead
        S_pq
        Y_n
    end
    DGUnitTrackingFRT(;
        K_pll,
        K_PT1,
        T_PT1,
        K_FRT,
        I_max,
        Vref,
        Vdead,
        S_pq,
        Y_n,
    ) = DGUnitTrackingFRT(
        K_pll,
        K_PT1,
        T_PT1,
        K_FRT,
        I_max,
        Vref,
        Vdead,
        S_pq,
        Y_n,
    )
    function construct_vertex(par::DGUnitTrackingFRT)
        K_pll = par.K_pll
        K_PT1 = par.K_PT1
        T_PT1 = par.T_PT1
        K_FRT = par.K_FRT
        I_max = par.I_max
        Vref = par.Vref
        Vdead = par.Vdead
        S_pq = par.S_pq
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            # current from incoming/outgoing edges, nodal shunt, PQ load background
            i = total_current(e_s, e_d) + Y_n * u - conj(S_pq) / conj(u)
            Vamp = abs(u)

            θ = x[3]
            P_int = x[4]
            Q_int = x[5]
            P_g = x[6]
            Q_g = x[7]

            # power flow tracking
            ΔP, ΔQ = p

            # detect error state
            over_voltage = Vamp > Vref + Vdead
            under_voltage = Vamp > Vref - Vdead
            smooth_fault_state(x) =
                SmoothStep(x, Vref + Vdead, Vref - Vdead; order = 100)

            # don't integrate the error during faults
            dP_int = -ΔP * smooth_fault_state(Vamp)
            dQ_int = -ΔQ * smooth_fault_state(Vamp)

            dP_g = (K_PT1 * P_int / 2.0 - P_g) / T_PT1
            dQ_g = (K_PT1 * Q_int / 2.0 - Q_g) / T_PT1

            # current set point adjustment
            I_r = complex(P_g, -Q_g) / Vamp

            # FRT
            # additional reactive current during fault
            iqplus = 0.0
            if over_voltage
                println("over")
                iqplus = -K_FRT * (Vamp - (Vref + Vdead))
            elseif under_voltage
                println("under")
                iqplus = -K_FRT * (Vamp - (Vref - Vdead))
            end

            I_r += complex(0.0, -iqplus)

            # PLL
            # (x[2] * cos(θ) - x[1] * sin(θ))
            dθ = -K_pll * Vamp * sin(θ - angle(u))

            du = i - I_r * exp(θ * 1im) # = 0 constraint

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                dx[3] = dθ
                dx[4] = dP_int
                dx[5] = dQ_int
                dx[6] = dP_g
                dx[7] = dQ_g
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(
            f! = rhs!,
            dim = 7,
            mass_matrix = Diagonal([0, 0, 1, 1, 1, 1, 1]),
            sym = Symbol[:u_r, :u_i, :θ, :P_int, :Q_int, :P_g, :Q_g],
        )
    end
    symbolsof(::DGUnitTrackingFRT) = begin
        [:u_r, :u_i, :θ, :P_int, :Q_int, :P_g, :Q_g]
    end
    dimension(::DGUnitTrackingFRT) = begin
        7
    end
end

begin
    @__doc__ struct StaticGenerator <: AbstractNode
        K_pll
        S_ref
        K_FRT
        I_max
        Vref
        Vdead
        Y_n
    end
    StaticGenerator(;
        K_pll,
        S_ref,
        K_FRT,
        I_max,
        Vref,
        Vdead,
        Y_n,
    ) = StaticGenerator(
        K_pll,
        S_ref,
        K_FRT,
        I_max,
        Vref,
        Vdead,
        Y_n,
    )
    function construct_vertex(par::StaticGenerator)
        K_pll = par.K_pll
        S_ref = par.S_ref
        K_FRT = par.K_FRT
        I_max = par.I_max
        Vref = par.Vref
        Vdead = par.Vdead
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            # current from incoming/outgoing edges, nodal shunt, PQ load background
            i = total_current(e_s, e_d) + Y_n(t) * u
            Vamp = abs(u)

            θ = x[3]
            i_d_r = x[4]
            i_q_r = x[5]

            # PLL
            dθ = -K_pll * Vamp * sin(θ - angle(u))

            # FRT
            # additional reactive current during fault
            if Vamp > Vref + Vdead
                iqplus = K_FRT * (Vamp - (Vref + Vdead)) # neg for over_voltage
            elseif Vamp < Vref - Vdead
                iqplus = K_FRT * (Vamp - (Vref - Vdead)) # pos for under_voltage
            else
                iqplus = 0.
            end

            # current set point adjustment
            I_r = conj(S_ref) / conj(u) #exp(θ * 1im) *

            # save values for plotting
            di_d_r = i_d_r - real(I_r)
            di_q_r = i_q_r - imag(I_r)

            # current limiter
            if Vamp > Vref + Vdead || Vamp < Vref - Vdead
                iq_fault = Limit(imag(I_r) + iqplus, -I_max, I_max)
                id_lim = sqrt(I_max^2 - iq_fault^2)
                id_fault = Limit(real(I_r), 0., id_lim) #-id_lim, id_lim)
                I_r_limited = complex(id_fault, iq_fault)
            else
                id_normal = Limit(real(I_r), 0., I_max) #-I_max, I_max)
                iq_lim = sqrt(I_max^2 - id_normal^2)
                iq_normal = Limit(imag(I_r), -iq_lim, iq_lim)
                I_r_limited = complex(id_normal, iq_normal)
            end

            du = i - I_r_limited # = 0 constraint

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                dx[3] = dθ
                dx[4] = di_d_r
                dx[5] = di_q_r
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(
            f! = rhs!,
            dim = 5,
            mass_matrix = Diagonal([0, 0, 1, 0, 0]),
            sym = Symbol[:u_r, :u_i, :θ, :i_d_r, :i_q_r],
        )
    end
    symbolsof(::StaticGenerator) = begin
        [:u_r, :u_i, :θ, :i_d_r, :i_q_r]
    end
    dimension(::StaticGenerator) = begin
        5
    end
end

function Limit(x, l, u)
    if x < l
        return l
    elseif x > u
        return u
    else
        return x
    end
end

begin
    @__doc__ struct StaticGeneratorObs <: AbstractNode
        K_pll
        S_ref
        K_FRT
        I_max
        Vref
        Vdead
        Y_n
    end
    StaticGeneratorObs(; K_pll, S_ref, K_FRT, I_max, Vref, Vdead, Y_n) =
        StaticGeneratorObs(K_pll, S_ref, K_FRT, I_max, Vref, Vdead, Y_n)
    function construct_vertex(par::StaticGeneratorObs)
        K_pll = par.K_pll
        S_ref = par.S_ref
        K_FRT = par.K_FRT
        I_max = par.I_max
        Vref = par.Vref
        Vdead = par.Vdead
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            # current from incoming/outgoing edges, nodal shunt, PQ load background
            i = total_current(e_s, e_d) + Y_n * u
            Vamp = abs(u)

            # this is the only dynamic variable
            θ = x[3]

            # to keep track of the internal currents for debugging ...
            i_d_r = x[4] # behind FRT
            i_q_r = x[5] # behind FRT
            i_d_r_b = x[6] # before FRT
            i_q_r_b = x[7] # before FRT
            iqp = x[8]

            # PLL
            dθ = -K_pll * Vamp * sin(θ - angle(u))

            # FRT
            # additional reactive current during fault
            if Vamp > Vref + Vdead
                iqplus = K_FRT * (Vamp - (Vref + Vdead)) # pos for over_voltage
            elseif Vamp < Vref - Vdead
                iqplus = K_FRT * (Vamp - (Vref - Vdead)) # neg for under_voltage
            else
                iqplus = 0.0
            end

            # current set point
            I_r = conj(S_ref) / Vamp

            # save values for plotting
            di_d_r_b = i_d_r_b - real(I_r)
            di_q_r_b = i_q_r_b - imag(I_r)
            diqp = iqp - iqplus

            # FRT limiter
            if Vamp > Vref + Vdead || Vamp < Vref - Vdead
                iq_fault = Limit(imag(I_r) + iqplus, -I_max, I_max)
                id_lim = sqrt(I_max^2 - iq_fault^2)
                id_fault = Limit(real(I_r), 0.0, id_lim) #-id_lim, id_lim)
                I_r_limited = complex(id_fault, iq_fault)
            else
                id_normal = Limit(real(I_r), 0.0, I_max) #-I_max, I_max)
                iq_lim = sqrt(I_max^2 - id_normal^2)
                iq_normal = Limit(imag(I_r), -iq_lim, iq_lim)
                I_r_limited = complex(id_normal, iq_normal)
            end

            # save values for plotting
            di_d_r = i_d_r - real(I_r_limited)
            di_q_r = i_q_r - imag(I_r_limited)

            du = i - exp(θ * 1im) * I_r_limited # = 0 constraint

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                dx[3] = dθ
                dx[4] = di_d_r
                dx[5] = di_q_r
                dx[6] = di_d_r_b
                dx[7] = di_q_r_b
                dx[8] = diqp
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(
            f! = rhs!,
            dim = 8,
            mass_matrix = Diagonal([0, 0, 1, 0, 0, 0, 0, 0]),
            sym = Symbol[
                :u_r,
                :u_i,
                :θ,
                :i_d_r,
                :i_q_r,
                :i_d_r_b,
                :i_q_r_b,
                :iqp,
            ],
        )
    end
    symbolsof(::StaticGeneratorObs) = begin
        [:u_r, :u_i, :θ, :i_d_r, :i_q_r, :i_d_r_b, :i_q_r_b, :iqp]
    end
    dimension(::StaticGeneratorObs) = begin
        8
    end
end

# x = -2:.01:2
# plot(x, SmoothStep.(x, 1, -1))
# plot!(x, SmoothClamp.(x, 1, -1))
function smooth_step(x, loc; order = 10)
    (1 .+ tanh.(order * (x - loc))) / 2.0
end

function smooth_cutoff(value, threshold; scale=10.)
    log(exp(scale * value) + exp.(scale * threshold)) / scale
end

begin
    @__doc__ struct DGUnit <: AbstractNode
        K_pll
        K_PT1
        T_PT1
        K_FRT
        I_max
        Vref
        Vdead
        S_pq
        Y_n
    end
    DGUnit(; K_pll, K_PT1, T_PT1, K_FRT, I_max, Vref, Vdead, S_pq, Y_n) =
        DGUnit(K_pll, K_PT1, T_PT1, K_FRT, I_max, Vref, Vdead, S_pq, Y_n)
    function construct_vertex(par::DGUnit)
        K_pll = par.K_pll
        K_PT1 = par.K_PT1
        T_PT1 = par.T_PT1
        K_FRT = par.K_FRT
        I_max = par.I_max
        Vref = par.Vref
        Vdead = par.Vdead
        S_pq = par.S_pq
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            # current from incoming/outgoing edges, nodal shunt, PQ load background
            i = total_current(e_s, e_d) + Y_n * u - conj(S_pq) / conj(u)
            Vamp = abs(u)

            # this are the dynamic variables
            θ = x[3]
            P_int = x[4]
            Q_int = x[5]
            P_g = x[6]
            Q_g = x[7]

            # to keep track of the internal currents for debugging ...
            iqp = x[8]

            # power flow tracking
            ΔP, ΔQ = p

            # detect error state
            over_voltage = Vamp > Vref + Vdead
            under_voltage = Vamp < Vref - Vdead
            # smooth_fault_state(x) = SmoothStep(x, Vref + Vdead, Vref - Vdead; order = 100)

            # don't integrate the error during faults
            if over_voltage | under_voltage
                dP_int = 0.
                dQ_int = 0.
                dP_g = 0.
                dQ_g = 0.
            else
                dP_int = -ΔP #* smooth_fault_state(Vamp)
                dQ_int = -ΔQ #* smooth_fault_state(Vamp)
                dP_g = (K_PT1 * P_int / 2.0 - P_g) / T_PT1
                dQ_g = (K_PT1 * Q_int / 2.0 - Q_g) / T_PT1
            end

            # dP_int = -ΔP
            # dQ_int = -ΔQ

            # dP_g = (K_PT1 * P_int / 2.0 - P_g) / T_PT1
            # dQ_g = (K_PT1 * Q_int / 2.0 - Q_g) / T_PT1

            # PLL
            dθ = -K_pll * Vamp * sin(θ - angle(u))

            # FRT
            # additional reactive current during fault
            if over_voltage
                iqplus = K_FRT * (Vamp - (Vref + Vdead)) # pos for over_voltage
            elseif under_voltage
                iqplus =-K_FRT * ((Vref - Vdead) - Vamp) # neg for under_voltage
            else
                iqplus = 0.0
            end

            # current set point
            I_r = conj(complex(P_g, Q_g)) / Vamp

            # save values for plotting
            diqp = iqp - iqplus

            # FRT limiter
            if over_voltage | under_voltage
                iq_fault = Limit(imag(I_r) + iqplus, -I_max, I_max)
                id_lim = sqrt(I_max^2 - iq_fault^2)
                id_fault = Limit(real(I_r), 0.0, id_lim) #-id_lim, id_lim)
                I_r_limited = complex(id_fault, iq_fault)
            else
                id_normal = Limit(real(I_r), 0.0, I_max) #-I_max, I_max)
                iq_lim = sqrt(I_max^2 - id_normal^2)
                iq_normal = Limit(imag(I_r), -iq_lim, iq_lim)
                I_r_limited = complex(id_normal, iq_normal)
            end

            du = i - exp(θ * 1im) * I_r_limited  # = 0 constraint

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                dx[3] = dθ
                dx[4] = dP_int
                dx[5] = dQ_int
                dx[6] = dP_g
                dx[7] = dQ_g
                dx[8] = diqp
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(
            f! = rhs!,
            dim = 8,
            mass_matrix = Diagonal([0, 0, 1, 1, 1, 1, 1, 0]),
            sym = Symbol[
                :u_r,
                :u_i,
                :θ,
                :P_int,
                :Q_int,
                :P_g,
                :Q_g,
                :iqp,
            ],
        )
    end
    symbolsof(::DGUnit) = begin
        [:u_r, :u_i, :θ, :P_int, :Q_int, :P_g, :Q_g, :iqp]
    end
    dimension(::DGUnit) = begin
        8
    end
end

function clamp(x, lowerlimit, upperlimit)
    if (x < lowerlimit)
        x = lowerlimit
    elseif (x > upperlimit)
        x = upperlimit
    end
    return x
end
function smootherstep(edge0, edge1, x)
    # Scale, and clamp x to 0..1 range
    x = clamp((x - edge0) / (edge1 - edge0), 0.0, 1.)
    # Evaluate polynomial
    return x * x * x * (x * (x * 6 - 15) + 10)
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
Threads.nthreads() ≥ 2 || error("At least two Julia threads required.")
println("--- :julia: Instantiating project")
using Pkg
Pkg.activate("..")
Pkg.instantiate()
Pkg.activate(".")
Pkg.instantiate()
push!(LOAD_PATH, joinpath(@__DIR__, ".."))
deleteat!(LOAD_PATH, 2)
println("+++ :julia: Building Literate.jl examples")

using ThreadPinning
using Literate

const src = "https://github.com/carstenbauer/ThreadPinning.jl"

cd(@__DIR__) do
    Literate.markdown(
        "src/examples/ex_core2core_latency.jl",
        "src/examples/";
        repo_root_url="$src/blob/main/docs",
        execute=true,
    )
end

ext: jl


text:
# This file was generated, do not modify it. # hide
y, X = unpack(hitters, ==(:Salary), col->true);


text:
module USTC_60

using Printf
using XEB
using OMEinsum, OMEinsumContractionOrders
using CUDA
using Graphs
using Statistics

vs_enable = setdiff(collect(0:65), [4, 5, 11, 17, 22, 65])
L1 = [0, 1, 2, 3, 6, 7, 8, 9, 10, 12, 13, 14, 15, 16, 18, 19, 20,
    21, 24, 25, 26, 30, 31, 32, 36, 37, 42, 43, 48, 54]
R1 = setdiff(vs_enable, L1)
L2 = [0, 1, 2, 3, 6, 7, 8, 9, 10, 12, 13, 14, 15, 18, 19, 20, 21, 24, 25,
    26, 27, 30, 31, 32, 33, 36, 37, 38, 42, 43, 44, 48, 49, 54, 55, 60]
R2 = setdiff(vs_enable, L2)
L3 = [0, 1, 2, 3, 6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 24, 25,
    26, 27, 30, 31, 32, 33, 36, 37, 38, 42, 43, 44, 48, 49, 54, 55, 60]
R3 = setdiff(vs_enable, L3)
L4 = [0, 1, 6, 7, 8, 12, 13, 14, 18, 19, 20, 21, 24, 25,
    26, 30, 31, 32, 36, 37, 42, 43, 48, 49, 54, 55, 60]
R4 = setdiff(vs_enable, L4)

L, R = L4, R4
D = 24
g = ustc_layout_60(60, D)
cuts = XEB.generate_cut(g, L, 1, D)
XEB.simplify!(g, cuts);
flux_g_L = 2.0^(flux(g, L) - length(L))
flux_g_R = 2.0^(flux(g, R) - length(R))
@show flux_g_L, flux_g_R
ec_L, ts_L, ids_size_L = to_ein_code(g, L; haar = true, enable_cuda = true)
ec_R, ts_R, ids_size_R = to_ein_code(g, R; haar = true, enable_cuda = true)
length(ec_L.iy), length(ec_R.iy)

ec_L_greedy = optimize_code(ec_L, ids_size_L, GreedyMethod())
ec_L_opt = optimize_code(ec_L_greedy, ids_size_L, TreeSA())
@show timespace_complexity(ec_L_opt, ids_size_L)
ec_R_greedy = optimize_code(ec_R, ids_size_R, GreedyMethod())
ec_R_opt = optimize_code(ec_R_greedy, ids_size_R, TreeSA())
@show timespace_complexity(ec_R_opt, ids_size_R)
CUDA.@time ec_L_opt(ts_L...);
CUDA.@time ec_R_opt(ts_R...);
xebs_L = Float64[]
xebs_R = Float64[]
maxs_L = Float64[]
maxs_R = Float64[]

for i = 1:100
    g = ustc_layout_60(60, D);
    cuts = XEB.generate_cut(g, L, 1, D)
    XEB.simplify!(g, cuts);
    _, ts_L, _ = to_ein_code(g, L; haar = false, enable_cuda = true)
    _, ts_R, _ = to_ein_code(g, R; haar = false, enable_cuda = true)
    p_L = ec_L_opt(ts_L...);
    p_R = ec_R_opt(ts_R...);
    # @show sum(p_L), sum(p_R)
    # p_L = p_L / sum(p_L);
    # p_R = p_R / sum(p_R);
    xeb_L = length(p_L)*sum(p_L .* p_L) - 1
    push!(xebs_L, xeb_L)
    xeb_R = length(p_R)*sum(p_R .* p_R) - 1
    push!(xebs_R, xeb_R)
    push!(maxs_L, length(p_L)*maximum(p_L) - 1)
    push!(maxs_R, length(p_R)*maximum(p_R) - 1)
    print("$(i): ")
    @printf "XEB_max = %.3e = %.3e + %.3e (flux = %.3e + %.3e)\n" (mean(maxs_L)+mean(maxs_R)) mean(maxs_L) mean(maxs_R) flux_g_L flux_g_R
    # println("$(i): XEB_max = $(mean(maxs_L)+mean(maxs_R)) = $(mean(maxs_L)) + $(mean(maxs_R))")
    # println("$(i): XEB_L = $(mean(xebs_L)), XEB_R = $(mean(xebs_R))")
end

end


text:
## A template for defining a new derivative evaluation method

# Define the method struct inheriting from either `GenerativeDerivativeMethod` or `NonGenerativeDerivativeMethod`
# Here we'll use `GenerativeDerivativeMethod` assuming our method will add extra supports
# This should contain any information that is needed to implement the method
struct MyDerivMethod <: GenerativeDerivativeMethod
    my_attr::Float64 # REPLACE WITH ACTUAL ATTRIBUTE
    # ADD ANY MORE INFORMATION THAT IS NEEDED
end

# Extend `generative_support_info` (only needed for generative methods)
function InfiniteOpt.generative_support_info(method::MyDerivMethod)
    info = UniformGenerativeInfo([method.my_attr], InternalLabel)
    return info # REPLACE WITH NEEDED AbstractGenerativeInfo THAT IS NEEDED TO MAKE THE GENERATIVE SUPPORTS
end

# Extend `evaluate_derivative`
# Generative methods must include calling `add_generative_supports`
# It will likely also be convenient to use `make_reduced_expr`
function InfiniteOpt.evaluate_derivative(
    dref::GeneralVariableRef,
    method::MyDerivMethod,
    write_model::JuMP.AbstractModel
    )::Vector{JuMP.AbstractJuMPScalar}
    # get the basic derivative information
    vref = derivative_argument(dref)
    pref = operator_parameter(dref)
    # make sure generative supports are added to the model
    InfiniteOpt.add_generative_supports(pref)
    # generate the derivative expressions h_i corresponding to equations of
    # the form h_i = 0 (REPLACE BELOW WITH ACTUAL IMPLEMENTATION)
    supps = supports(pref, label = All)
    exprs = Vector{JuMP.AbstractJuMPScalar}(undef, length(supps) - 1)
    for i in eachindex(exprs)
        d = InfiniteOpt.make_reduced_expr(dref, pref, supps[i], write_model)
        v1 = InfiniteOpt.make_reduced_expr(vref, pref, supps[i], write_model)
        v2 = InfiniteOpt.make_reduced_expr(vref, pref, supps[i + 1], write_model)
        change = supps[i + 1] - supps[i]
        exprs[i] = JuMP.@expression(write_model, -change * d + v2 - v1)
    end
    return exprs
end


text:
import Base: rand,
             ∈

export Ballp

"""
    Ballp{N<:AbstractFloat, VN<:AbstractVector{N}} <: AbstractCentrallySymmetric{N}

Type that represents a ball in the p-norm, for ``1 ≤ p ≤ ∞``.

It is defined as the set

```math
\\mathcal{B}_p^n(c, r) = \\{ x ∈ \\mathbb{R}^n : ‖ x - c ‖_p ≤ r \\},
```
where ``c ∈ \\mathbb{R}^n`` is its center and ``r ∈ \\mathbb{R}_+`` its radius.
Here ``‖ ⋅ ‖_p`` for ``1 ≤ p ≤ ∞`` denotes the vector ``p``-norm, defined as
``‖ x ‖_p = \\left( \\sum\\limits_{i=1}^n |x_i|^p \\right)^{1/p}`` for any
``x ∈ \\mathbb{R}^n``.

### Fields

- `p`      -- norm as a real scalar
- `center` -- center of the ball as a real vector
- `radius` -- radius of the ball as a scalar (``≥ 0``)

## Notes

The special cases ``p=1``, ``p=2`` and ``p=∞`` fall back to the specialized
types `Ball1`, `Ball2` and `BallInf`, respectively.

### Examples

A five-dimensional ball in the ``p=3/2`` norm centered at the origin of radius
0.5:

```jldoctest ballp_constructor
julia> B = Ballp(3/2, zeros(5), 0.5)
Ballp{Float64,Array{Float64,1}}(1.5, [0.0, 0.0, 0.0, 0.0, 0.0], 0.5)

julia> dim(B)
5
```

We evaluate the support vector in direction ``[1,2,…,5]``:

```jldoctest ballp_constructor
julia> σ([1., 2, 3, 4, 5], B)
5-element Array{Float64,1}:
 0.013516004434607558
 0.05406401773843023
 0.12164403991146802
 0.21625607095372093
 0.33790011086518895
```
"""
struct Ballp{N<:AbstractFloat, VN<:AbstractVector{N}} <: AbstractCentrallySymmetric{N}
    p::N
    center::VN
    radius::N

    # default constructor with domain constraint for radius and p
    function Ballp(p::N, center::VN, radius::N) where {N<:Real, VN<:AbstractVector{N}}
        @assert radius >= zero(N) "radius must not be negative"
        @assert p >= 1 "p must not be less than 1"
        if p == Inf
            return BallInf(center, radius)
        elseif p == 2
            return Ball2(center, radius)
        elseif p == 1
            return Ball1(center, radius)
        else
            return new{N, VN}(p, center, radius)
        end
    end
end

isoperationtype(::Type{<:Ballp}) = false
isconvextype(::Type{<:Ballp}) = true

# --- AbstractCentrallySymmetric interface functions ---

"""
    center(B::Ballp{N}) where {N<:AbstractFloat}

Return the center of a ball in the p-norm.

### Input

- `B` -- ball in the p-norm

### Output

The center of the ball in the p-norm.
"""
function center(B::Ballp{N}) where {N<:AbstractFloat}
    return B.center
end

# --- LazySet interface functions ---

"""
    σ(d::AbstractVector{N}, B::Ballp{N}) where {N<:AbstractFloat}

Return the support vector of a `Ballp` in a given direction.

### Input

- `d` -- direction
- `B` -- ball in the p-norm

### Output

The support vector in the given direction.
If the direction has norm zero, the center of the ball is returned.

### Algorithm

The support vector of the unit ball in the ``p``-norm along direction ``d`` is:

```math
σ_{\\mathcal{B}_p^n(0, 1)}(d) = \\dfrac{\\tilde{v}}{‖\\tilde{v}‖_q},
```
where ``\\tilde{v}_i = \\frac{|d_i|^q}{d_i}`` if ``d_i ≠ 0`` and
``tilde{v}_i = 0`` otherwise, for all ``i=1,…,n``, and ``q`` is the conjugate
number of ``p``.
By the affine transformation ``x = r\\tilde{x} + c``, one obtains that
the support vector of ``\\mathcal{B}_p^n(c, r)`` is

```math
σ_{\\mathcal{B}_p^n(c, r)}(d) = \\dfrac{v}{‖v‖_q},
```
where ``v_i = c_i + r\\frac{|d_i|^q}{d_i}`` if ``d_i ≠ 0`` and ``v_i = 0``
otherwise, for all ``i = 1, …, n``.
"""
function σ(d::AbstractVector{N}, B::Ballp{N}) where {N<:AbstractFloat}
    p = B.p
    q = p/(p-1)
    v = similar(d)
    @inbounds for (i, di) in enumerate(d)
        v[i] = di == zero(N) ? di : abs.(di).^q / di
    end
    vnorm = norm(v, p)
    svec = vnorm != zero(N) ? @.(B.center + B.radius * (v/vnorm)) : B.center
    return svec
end

"""
    ∈(x::AbstractVector{N}, B::Ballp{N}) where {N<:AbstractFloat}

Check whether a given point is contained in a ball in the p-norm.

### Input

- `x` -- point/vector
- `B` -- ball in the p-norm

### Output

`true` iff ``x ∈ B``.

### Notes

This implementation is worst-case optimized, i.e., it is optimistic and first
computes (see below) the whole sum before comparing to the radius.
In applications where the point is typically far away from the ball, a fail-fast
implementation with interleaved comparisons could be more efficient.

### Algorithm

Let ``B`` be an ``n``-dimensional ball in the p-norm with radius ``r`` and let
``c_i`` and ``x_i`` be the ball's center and the vector ``x`` in dimension
``i``, respectively.
Then ``x ∈ B`` iff ``\\left( ∑_{i=1}^n |c_i - x_i|^p \\right)^{1/p} ≤ r``.

### Examples

```jldoctest
julia> B = Ballp(1.5, [1., 1.], 1.)
Ballp{Float64,Array{Float64,1}}(1.5, [1.0, 1.0], 1.0)

julia> [.5, -.5] ∈ B
false

julia> [.5, 1.5] ∈ B
true
```
"""
function ∈(x::AbstractVector{N}, B::Ballp{N}) where {N<:AbstractFloat}
    @assert length(x) == dim(B)
    sum = zero(N)
    for i in eachindex(x)
        sum += abs(B.center[i] - x[i])^B.p
    end
    return sum^(one(N)/B.p) <= B.radius
end

"""
    rand(::Type{Ballp}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)

Create a random ball in the p-norm.

### Input

- `Ballp` -- type for dispatch
- `N`     -- (optional, default: `Float64`) numeric type
- `dim`   -- (optional, default: 2) dimension
- `rng`   -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`  -- (optional, default: `nothing`) seed for reseeding

### Output

A random ball in the p-norm.

### Algorithm

The center and radius are normally distributed with mean 0 and standard
deviation 1.
Additionally, the radius is nonnegative.
The p-norm is a normally distributed number ≥ 1 with mean 1 and standard
deviation 1.
"""
function rand(::Type{Ballp};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing)
    rng = reseed(rng, seed)
    p = one(N) + abs(randn(rng, N))
    center = randn(rng, N, dim)
    radius = abs(randn(rng, N))
    return Ballp(p, center, radius)
end

"""
    translate(B::Ballp{N}, v::AbstractVector{N}) where {N<:AbstractFloat}

Translate (i.e., shift) a ball in the p-norm by a given vector.

### Input

- `B` -- ball in the p-norm
- `v` -- translation vector

### Output

A translated ball in the p- norm.

### Algorithm

We add the vector to the center of the ball.
"""
function translate(B::Ballp{N}, v::AbstractVector{N}) where {N<:AbstractFloat}
    @assert length(v) == dim(B) "cannot translate a $(dim(B))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    return Ballp(B.p, center(B) + v, B.radius)
end


text:
function generate_abc_sample_serial(plan::abc_pmc_plan_type, sampler::Distribution, ss_true, epsilon::Real; n::Integer = 100)
  #num_param = length(Distributions.rand(plan.prior))
  num_param = length(plan.prior)
  theta_plot = Array{Float64}(undef,num_param, n)
  dist_plot = Array{Float64}(undef,n)
  for i in 1:n
    theta_plot[:,i], dist_plot[i], attempts_plot, accept_log_plot, reject_log_plot = ABC.generate_theta(abc_plan, sampler, ss_true, epsilon)
  end
  return theta_plot, dist_plot
end

# Generate initial abc population from prior, aiming for d(ss,ss_true)<epsilon
function init_abc_serial(plan::abc_pmc_plan_type, ss_true)
  num_param = length(Distributions.rand(plan.prior))
  # Allocate arrays
  theta = Array{Float64}(undef,num_param,plan.num_part)
  dist_theta = Array{Float64}(undef,plan.num_part)
  attempts = zeros(plan.num_part)
  accept_log_combo = abc_log_type()
  reject_log_combo = abc_log_type()
  # Draw initial set of theta's from prior (either with dist<epsilon or best subject to num_max_attempts and num_max_valid_attempts)

  if plan.save_params || plan.save_summary_stats || plan.save_distances
     push!(accept_log_combo.generation_starts_at,length(accept_log_combo.dist)+1)
     push!(reject_log_combo.generation_starts_at,length(reject_log_combo.dist)+1)
  end
  for i in 1:plan.num_part
      theta[:,i], dist_theta[i], attempts[i], accept_log, reject_log = generate_theta(plan, plan.prior, ss_true, plan.epsilon_init)
      append_to_abc_log!(accept_log_combo,plan,accept_log.theta,accept_log.ss,accept_log.dist)
      append_to_abc_log!(reject_log_combo,plan,reject_log.theta,reject_log.ss,reject_log.dist)
  end
  weights = fill(1.0/plan.num_part,plan.num_part)
  logpriorvals = Distributions.logpdf(plan.prior,theta)
  push!(accept_log_combo.epsilon,plan.epsilon_init)
  push!(reject_log_combo.epsilon,plan.epsilon_init)

  return abc_population_type(theta,weights,dist_theta,maximum(dist_theta),logpriorvals,accept_log_combo,reject_log_combo)
end

# Update the abc population once
function update_abc_pop_serial(plan::abc_pmc_plan_type, ss_true, pop::abc_population_type, sampler::Distribution, epsilon::Float64;
                        attempts::Array{Int64,1} = zeros(Int64,plan.num_part) )
  #new_pop = deepcopy(pop)
  new_pop = abc_population_type(size(pop.theta,1), size(pop.theta,2), accept_log=pop.accept_log, reject_log=pop.reject_log, repeats=pop.repeats)
  if plan.save_params || plan.save_summary_stats || plan.save_distances
     push!(new_pop.accept_log.generation_starts_at,length(new_pop.accept_log.dist)+1)
     push!(new_pop.reject_log.generation_starts_at,length(new_pop.reject_log.dist)+1)
  end
  median_logpdf_old = median(pop.logpdf)
     for i in 1:plan.num_part
       theta_star, dist_theta_star, attempts[i], accept_log, reject_log = generate_theta(plan, sampler, ss_true, epsilon)
       append_to_abc_log!(new_pop.accept_log,plan,accept_log.theta,accept_log.ss,accept_log.dist)
       append_to_abc_log!(new_pop.reject_log,plan,reject_log.theta,reject_log.ss,reject_log.dist)

       #if dist_theta_star < epsilon # replace theta with new set of parameters and update weight
         @inbounds new_pop.dist[i] = dist_theta_star
         @inbounds new_pop.theta[:,i] = theta_star
         prior_logpdf = Distributions.logpdf(plan.prior,theta_star)
         #if isa(prior_logpdf, Array)   # TODO: Can remove this once Danley's code uses composite distribution
         #   prior_logpdf = prior_logpdf[1]
         #end
         # sampler_pdf calculation must match distribution used to update particle
         sampler_logpdf = logpdf(sampler, theta_star )
         #@inbounds new_pop.weights[i] = prior_pdf/sampler_pdf
         @inbounds new_pop.weights[i] = exp(prior_logpdf-sampler_logpdf)
         @inbounds new_pop.logpdf[i] = sampler_logpdf
      if dist_theta_star < epsilon # replace theta with new set of parameters and update weight
         @inbounds new_pop.repeats[i] = 0
       else  # failed to generate a closer set of parameters, so...
         # ... generate new data set with existing parameters
         #new_data = plan.gen_data(theta_star)
         #new_ss = plan.calc_summary_stats(new_data)
         #@inbounds new_pop.dist[i] = plan.calc_dist(ss_true,new_ss)
         # ... just keep last value for this time, and mark it as a repeat
         #=
         @inbounds new_pop.dist[i] = pop.dist[i]
         @inbounds new_pop.theta[:,i] = pop.theta[:,i]
         @inbounds new_pop.weights[i] = pop.weights[i]
         @inbounds new_pop.logpdf[i] = pop.logpdf[i]
         =#
         @inbounds new_pop.repeats[i] += 1
       end
     end # i / num_parts
   new_pop.weights ./= sum(new_pop.weights)
   push!(new_pop.accept_log.epsilon,epsilon)
   push!(new_pop.reject_log.epsilon,epsilon)
   #println("New pop weights = ", new_pop.weights)
   return new_pop
end

#=
# run the ABC algorithm matching to summary statistics ss_true, starting from an initial population (e.g., output of previosu call)
function run_abc(plan::abc_pmc_plan_type, ss_true, pop::abc_population_type; verbose::Bool = false, print_every::Integer=1, in_parallel::Bool = plan.in_parallel )
  attempts = zeros(Int64,plan.num_part)
  # Set initial epsilon tolerance based on current population
  epsilon = quantile(pop.dist,plan.init_epsilon_quantile)
  eps_diff_count = 0
  for t in 1:plan.num_max_times
    local new_pop
    sampler = plan.make_proposal_dist(pop, plan.tau_factor)
    if plan.adaptive_quantiles
      min_quantile = t==1 ? 4.0/size(pop.theta,2) : 0.5 # 1.0/sqrt(size(pop.theta,2))
      epsilon = choose_epsilon_adaptive(pop, sampler, min_quantile=min_quantile)
    end
    if in_parallel
      #new_pop = update_abc_pop_parallel_darray(plan, ss_true, pop, epsilon, attempts=attempts)
       new_pop = update_abc_pop_parallel_distributed_map(plan, ss_true, pop, sampler, epsilon, attempts=attempts)
    else
	  #new_pop = update_abc_pop_serial(plan, ss_true, pop, epsilon, attempts=attempts)
      new_pop = update_abc_pop_serial(plan, ss_true, pop, sampler, epsilon, attempts=attempts)
    end
    pop = new_pop
    if verbose && (t%print_every == 0)
       println("# t= ",t, " eps= ",epsilon, " med(d)= ",median(pop.dist), " attempts= ",median(attempts), " ",maximum(attempts), " reps= ", sum(pop.repeats), " ess= ",ess(pop.weights)) #," mean(theta)= ",mean(pop.theta,2) )#) #, " tau= ",diag(tau) ) #
       # println("# t= ",t, " eps= ",epsilon, " med(d)= ",median(pop.dist), " max(d)= ", maximum(pop.dist), " med(attempts)= ",median(attempts), " max(a)= ",maximum(attempts), " reps= ", sum(pop.repeats), " ess= ",ess(pop.weights,pop.repeats)) #," mean(theta)= ",mean(pop.theta,2) )#) #, " tau= ",diag(tau) ) #
       println("# Mean(theta)= ", mean(pop.theta, 2), " Stand. Dev.(theta)= ", std(pop.theta, 2))
    end
    #if epsilon < plan.target_epsilon  # stop once acheive goal
    if maximum(pop.dist) < plan.target_epsilon  # stop once acheive goal
       println("# Reached ",epsilon," after ", t, " generations.")
       break
    end
    if median(attempts)>0.2*plan.num_max_attempt
      println("# Halting due to ", median(attempts), " median number of valid attempts.")
      break
    end
    if sum(pop.repeats)>plan.num_part
      println("# Halting due to ", sum(pop.repeats), " repeats. (Due to change, this is poor word choice. Really, it's too many failed attempts to generate parameters resulting in distance less than epsilon.)")
      break
    end
    eps_old = epsilon
    if maximum(attempts)<0.75*plan.num_max_attempt
      epsilon = minimum([maximum(pop.dist),epsilon * plan.epsilon_reduction_factor])
    end
    if ((abs(eps_old-epsilon)/epsilon) < 1.0e-5)
      eps_diff_count += 1
    else
      eps_diff_count = 0
    end
    if eps_diff_count > 1
      println("# Halting due to epsilon not improving significantly for 3 consecutive generations.")
      break
    end
  end # t / num_times
 #println("mean(theta) = ",[ sum(pop.theta[i,:])/size(pop.theta,2) for i in 1:size(pop.theta,1) ])
  if verbose
     println("# Epsilon history = ", pop.accept_log.epsilon)
     #println("Epsilon history = ", eps_arr)
     #println("Mean history = ", mean_arr)
     #println("Std Dev. history = ", std_arr)
  end
  return pop
end
=#


text:
#-----------------------------------------------------------------#
# This routine advances the solution in time using Midpoint RK2
# Written by F.X. Giraldo / P.R. Mugg on 9/19/2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey, CA 93943-5216
#------------------------------------------------------------------#

function rk2(q0,c,Δt,time)

    qi = q0

    #Evaluate at Beginning of Interval
    (rhs,f,g) = rhs_function(q0,c,time)

    #update to Midpoint of Interval
    qi = q0 + (Δt/2)*rhs

    #update to End of Interval
    (rhs,f,g) = rhs_function(qi,c,time+(Δt/2))
    qi = q0 + Δt*rhs

    #update
    return qp = qi
end


text:
function e001(N)
    return sum_threes_and_fives(N)
end

function sum_threes_and_fives(N)
    s = 0
    for i = 3:3:N-1
        s += i
    end
    for i = 5:5:N-1
        if i % 3 != 0
            s += i
        end
    end
    return s
end

function solution001()
    N = 1000
    answer001 = e001(N)
    println("The answer to problem 001 is $answer001")
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
using Gadfly, DataArrays, RDatasets, Distributions, Base.Test

set_default_plot_size(6inch, 3inch)

global density_dark_tested

p = Gadfly.with_theme(:dark) do
    plot(dataset("ggplot2", "diamonds"), x="Price", color="Cut", Geom.density)
end

# prevent these tests from running more than once
if !isdefined(:density_dark_tested)
    svg_str_dark = stringmime(MIME("image/svg+xml"), p)
    @test contains(svg_str_dark, hex(Gadfly.dark_theme.default_color))
    @test contains(svg_str_dark, hex(Gadfly.dark_theme.background_color))
    @test contains(svg_str_dark, hex(Gadfly.dark_theme.panel_fill))

    # Test reset.
    p2 = plot(dataset("ggplot2", "diamonds"), x="Price", color="Cut", Geom.density)
    svg_str_light = stringmime(MIME("image/svg+xml"), p2)
    @test !contains(svg_str_light, hex(Gadfly.dark_theme.default_color))
    @test !contains(svg_str_light, hex(Gadfly.dark_theme.background_color))
    @test !contains(svg_str_light, hex(Gadfly.dark_theme.panel_fill))

    density_dark_tested=true
end

p


text:
@inline function skycolor(ray::Ray{T}) where T
	white = SA[1.0, 1.0, 1.0]
	skyblue = SA[0.5, 0.7, 1.0]
	t = T(0.5)*(ray.dir.y + one(T))
    (one(T)-t)*white + t*skyblue
end

@inline color_vec3_in_rgb(v::Vec3{T}) where T = 0.5normalize(v) + SA{T}[0.5,0.5,0.5]

"""Compute color for a ray, recursively

	Args:
		depth: how many more levels of recursive ray bounces can we still compute?"""
@inline @fastmath function ray_color(r::Ray{T}, world::HittableList, depth=16) where T
    if depth <= 0
		return SA{T}[0,0,0]
	end

	rec::Union{HitRecord,Nothing} = hit(world, r, T(1e-4), typemax(T))
    if rec !== nothing
		# For debugging, represent vectors as RGB:
		# claforte TODO: adapt to latest code!
		# return color_vec3_in_rgb(rec.p) # show the normalized hit point
		# return color_vec3_in_rgb(rec.n⃗) # show the normal in RGB
		# return color_vec3_in_rgb(rec.p + rec.n⃗)
		# return color_vec3_in_rgb(random_vec3_in_sphere())
		#return color_vec3_in_rgb(rec.n⃗ + random_vec3_in_sphere())

        s = scatter(rec.mat, r, rec)
		if s.reflected
			return s.attenuation .* ray_color(s.r, world, depth-1)
		else
			return SA{T}[0,0,0]
		end
    else
        skycolor(r)
    end
end


text:
"""
    calmax(x)

Specifies maximum element for display purposes. Defaults to the maximum of `x`.
"""
calmax(x::Any) = getter(x, :calmax, i -> _caltype(x), i -> _calmax(i))
calmax!(x::Any, val::Any) = setter!(x, :calmax, val, i -> _caltype(i))

"""
    calmin(x)

Specifies minimum element for display purposes. Defaults to the minimum of `x`.
"""
calmin(x) = getter(x, :calmin, i -> _caltype(i), i -> _calmin(i))
calmin!(x, val) = setter!(x, :calmin, val, i -> _caltype(i))

###
_caltype(x::AbstractArray{T}) where {T} = T
_caltype(x::Any) = Float64
_calmax(x::AbstractArray) = maximum(x)
_calmax(x::Any) = one(Float64)
_calmin(x::AbstractArray) = minimum(x)
_calmin(x::Any) = one(Float64)

"""
    description(x) -> String

Retrieves description field that may say whatever you like.
"""
description(x) = getter(x, "description", String, i -> "")
description!(x, val) = setter!(x, "description", String, val)

"""
    freqdim(x) -> Int

Which spatial dimension (1, 2, or 3) corresponds to phase acquisition. If not
applicable to scan type defaults to `1`.
"""
freqdim(x) = getter(x, :freqdim, Int, 1)
freqdim!(x, val) = setter!(x, :freqdim, Int, val)

"""
    is_anatomical(::T) -> Bool

Returns `true` if `T` represents anatomical data.
"""
is_anatomical(::T) where {T} = is_anatomical(T)
is_anatomical(::Type{T}) where {T} = false

"""
    is_electrophysiology(::T) -> Bool

Returns `true` if `T` represents electrophysiology data.
"""
is_electrophysiology(::T) where {T} = is_electrophysiology(T)
is_electrophysiology(::Type{T}) where {T} = false

"""
    is_functional(::T) -> Bool

Returns `true` if `T` represents functional data.
"""
is_functional(::T) where {T} = is_functional(T)
is_functional(::Type{T}) where {T} = false
"""
    phasedim(x) -> Int

Which spatial dimension (1, 2, or 3) corresponds to phase acquisition.
"""
phasedim(x) = getter(x, :phasedim, Int, i -> 1)
phasedim!(x, val) = setter!(x, :phasedim, Int, val)

"""
    slicedim(x) -> Int

Which dimension slices where acquired at throughout MRI acquisition.
"""
slicedim(x) = getter(x, :slicedim, Int, i -> 1)
slicedim!(x, val) = setter!(x, :slicedim, Int, val)

"""
    slice_start(x) -> Int

Which slice corresponds to the first slice acquired during MRI acquisition
(i.e. not padded slices). Defaults to `1`.
"""
slice_start(x) = getter(x, :slice_start, Int, i -> 1)
slice_start!(x, val) = setter!(x, :slice_start, Int, val)

"""
    slice_end(x) -> Int

Which slice corresponds to the last slice acquired during MRI acquisition
(i.e. not padded slices).
"""
slice_end(x) = getter(x, :slice_end, Int, i -> 1)
slice_end!(x, val) = setter!(x, :slice_end, Int, val)

### TODO: should these go in ImageCore.jl ?
"""
    spatial_offset(x)

Provides the offset of each dimension (i.e., where each spatial axis starts).
"""
spatial_offset(x) = first.(coords_indices(x))

"""
    spatial_units(x)

Returns the units (i.e. Unitful.unit) that each spatial axis is measured in. If not
available `nothing` is returned for each spatial axis.
"""
spatial_units(x) = unit.(eltype.(spatial_indices(x)))

"""
    time_units(x)

Returns the units (i.e. Unitful.unit) the time axis is measured in. If not available
`nothing` is returned.
"""
time_units(x) = unit(eltype(timeaxis(x)))

"""
    stop_time(x)

Returns start time in seconds in relation to the stop of acquisition of the first
data sample in the corresponding neural dataset.
"""
stop_time(x) = last(timeaxis(x))

""""
    duration(x)

Duration of the event along the time axis in seconds.
"""
duration(x) = stop_time(x) - start_time(x)


text:
# Autogenerated wrapper script for H3_jll for x86_64-unknown-freebsd
export libh3

JLLWrappers.@generate_wrapper_header("H3")
JLLWrappers.@declare_library_product(libh3, "libh3.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libh3,
        "lib/libh3.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "fragments" begin

D = ParserCombinator.Parsers.DOT

for s in ("", " ", "  ", " // ", " /*  */ ", "\n\t")
    parse_one(s, Trace(D.spc_star + Eos()))
end

# test wrd too for high bit

str_one = D.str_one > D.unesc_join
@test parse_one("\"abc\"", str_one)[1] == "abc"
@test parse_one("\"abc\\\ndef\"", str_one)[1] == "abcdef"
@test parse_one("\"abc\\\ndef\\\nghi\"", str_one)[1] == "abcdefghi"
@test parse_one("\"abc\\\n\"", str_one)[1] == "abc"
@test parse_one("\"a\\\"c\\\n\"", str_one)[1] == "a\"c"
@test parse_one("\"\\\ndef\"", str_one)[1] == "def"
@test parse_one("\"\\\n\"", str_one)[1] == ""

@test parse_one("\"a\" +\"b\\\nc\"", D.str_many)[1] == "abc"

@test parse_one("\"a\" +\"b\\\nc\"", D.str_id)[1] == StringID("abc")
@test parse_one("abc", D.str_id)[1] == StringID("abc")

@test parse_one("-3.14", D.num_id)[1] == NumericID("-3.14")

@test parse_one("<<abc/>>", D.html_id)[1] == HtmlID("<abc/>")
@test parse_one("<<abc></abc>>", D.html_id)[1] == HtmlID("<abc></abc>")
@test parse_one("<<abc> \n</abc>>", D.html_id)[1] == HtmlID("<abc> \n</abc>")
@test parse_one("<<abc><def/></abc>>", D.html_id)[1] == HtmlID("<abc><def/></abc>")
@test parse_one("<<abc><abc></abc></abc>>", D.html_id)[1] == HtmlID("<abc><abc></abc></abc>")

@test parse_one("abc", D.id)[1] == StringID("abc")
@test parse_one("\"abc\"", D.id)[1] == StringID("abc")
@test parse_one("123", D.id)[1] == NumericID("123")
@test parse_one("<<abc/>>", D.id)[1] == HtmlID("<abc/>")

@test parse_one(":A:n", D.port)[1] == Port(StringID("A"), "n")
@test parse_one(":n", D.port)[1] == Port("n")
@test parse_one(":A", D.port)[1] == Port(StringID("A"))

@test parse_one("[a=b]", D.attr_list)[1] ==
Attribute[Attribute(StringID("a"), StringID("b"))]
@test parse_one("[a=b c = d]", D.attr_list)[1] ==
Attribute[Attribute(StringID("a"), StringID("b")),
          Attribute(StringID("c"), StringID("d"))]
@test parse_one("[a=b; c = d]", D.attr_list)[1] ==
Attribute[Attribute(StringID("a"), StringID("b")),
          Attribute(StringID("c"), StringID("d"))]
@test parse_one("[a=b,c=d][e=f]", D.attr_list)[1] ==
Attribute[Attribute(StringID("a"), StringID("b")),
          Attribute(StringID("c"), StringID("d")),
          Attribute(StringID("e"), StringID("f"))]

@test parse_one("a:b:c[d=e]", D.node_stmt)[1] ==
Node(NodeID(StringID("a"), Port(StringID("b"), "c")),
     Attribute[Attribute(StringID("d"), StringID("e"))])

@test parse_one("a--b[c=d]", D.edge_stmt)[1] ==
Edge(EdgeNode[NodeID(StringID("a")), NodeID(StringID("b"))],
     Attribute[Attribute(StringID("c"), StringID("d"))])

@test parse_one("graph [a=b]", D.attr_stmt)[1] ==
GraphAttributes(Attribute[Attribute(StringID("a"), StringID("b"))])

println("fragments ok")

end


text:
using MiniFB

function onclick(window, button, mod, isPressed)::Cvoid
    if Bool(isPressed)
        println("mouse clicked")
    end
    mousex = mfb_get_mouse_x(window)
    mousey = mfb_get_mouse_y(window)
    println("x: $mousex y: $mousey")
end

function active_fn(win::Ptr{Cvoid}, is_active::Bool)
    println("Window is now ",  is_active ? "" : "in", "active")
end

macro play(w, h, body)
    quote
        window = mfb_open_ex("Luxor -> Julia", $w, $h, MiniFB.WF_RESIZABLE)

        mfb_set_active_callback(window, active_fn)
        mfb_set_mouse_button_callback(window, onclick)

        buffer = zeros(UInt32, $w, $h)
        while true
            Drawing($w, $h, :image)
            origin()
            $(esc(body))
            m = permutedims(image_as_matrix!(buffer), (2, 1))
            finish()
            state = mfb_update(window, m)
            if state != MiniFB.STATE_OK
                break
            end
        end
        mfb_close(window)
    end
end

#= Examples:

1 clock

using Luxor, Colors, Dates, ColorSchemes

include(dirname(pathof(Luxor)) * "/play.jl")

function clock(cscheme=ColorSchemes.leonardo)
    @play 400 600 begin
        fontface("JuliaMono-Regular")

        sethue(get(cscheme, .0))
        paint()
        fontsize(30)

        sethue(get(cscheme, .2))
        h = Dates.hour(now())
        sector(O, 180, 200, 3π/2, 3π/2 + rescale(h, 0, 24, 0, 2pi), :fill)

        sethue(get(cscheme, .4))
        m = Dates.minute(now())
        sector(O, 160, 180, 3π/2, 3π/2 + rescale(m, 0, 60, 0, 2pi), :fill)

        sethue(get(cscheme, .6))
        s = Dates.second(now())
        sector(O, 140, 160, 3π/2, 3π/2 + rescale(s, 0, 60, 0, 2pi), :fill)

        sethue(get(cscheme, .8))
        ms = Dates.value(Dates.Millisecond(Dates.now()))
        sector(O, 137, 140, 3π/2, 3π/2 + rescale(ms, 0, 1000, 0, 2pi), :fill)

        sethue(get(cscheme, 1.0))
        text(Dates.format(Dates.now(), "HH:MM:SS"), halign=:center)
    end
end

clock(ColorSchemes.botticelli)

=#

#=

2: bezier path thingy

using Luxor, Colors, Dates

function bez()
    j = rand(1:15)
    n = 1
    pgon = [Point(-100, -100), Point(100, 100), Point(100, -50), Point(-50, 100)]
    @play 500 500 begin
        t1 = time_ns()
        background(0, 0, 0)
        sethue(HSB(285, .7, .9))
        k = 3
        #        for i in 1:5
        for s in 1:30
            setopacity(rescale(s, 1, 30, 0.1, 1.0))
            pgon[1] += (k * rand(-1:1), k * rand(-1:1))
            pgon[2] += (k * rand(-1:1), k * rand(-1:1))
            pgon[3] += (k * rand(-1:1), k * rand(-1:1))
            pgon[4] += (k * rand(-1:1), k * rand(-1:1))
            bp = makebezierpath(pgon)
            drawbezierpath(bp, :stroke)
        end
        #        end
        if !isinside(polycentroid(pgon), box(BoundingBox() * 0.8), allowonedge=true)
            pgon = [Point(-100, -100), Point(100, 100), Point(100, -50), Point(-50, 100)]
        end

        n = mod1(n + 1, 200)

        sethue("white")
        fontface("JuliaMono-Black")
        fontsize(20)
        text(Dates.format(now(), "HH:MM:SS"), halign=:center)

        # time calc
        sethue("white")
        t2 = time_ns()
        text(string("fps: ", round(1E9/(t2-t1), digits=2)), boxtopleft(BoundingBox() * 0.9))
        t1 = t2
        sleep(0.1)
    end
end

bez()

=#

#=

3 some balls

using Luxor, Colors, Dates

mutable struct Ball
    position::Point
    velocity::Point
    radius::Float64
    hue::Float64
end

function f()
    balls = [Ball(rand(BoundingBox(Point(-300, -300), Point(300, 300))), rand(BoundingBox(Point(-300, -300), Point(300, 300))), rand(1:10), rand(0:360)) for i in 1:450]
    @play 600 600 begin
        w = 500
        h = 500
        background("black")

        for ball in balls
            sethue(HSB(ball.hue, .8, 0.8))
            circle(ball.position, ball.radius, :fill)

            if !(-w/2 < ball.position.x < w/2)
                ball.velocity = -ball.velocity
                ball.position = Point(-w/2, ball.position.y)
            end
            if !(-h/2 < ball.position.y < h/2)
                ball.velocity = -ball.velocity
                ball.position = Point(ball.position.x, -h/2)
            end
            ball.velocity = Point(rand(-4:4, 2) ...)

            ball.position = ball.position + ball.velocity
        end
        sleep(0.01)
    end
end

f()
=#


text:
function add_data!(v::VegaVisualization;
                              name::AbstractString = string("table_", lowercase(randstring(5))),
                              x::AbstractVector = Int[],
                              x2::AbstractVector = Int[],
                              y::AbstractVector = Int[],
                              y2::AbstractVector = Int[],
                              group::AbstractVector = Int[])

    #Check arrays all of same length
    n = length(x)

    if n != length(y) || (n != length(group) && length(group) != 0)
        error("x, y and group must all have the same length")
    end

    #Create empty arrays as necessary
    res = Array(Dict{Any, Any}, n)
    isempty(group)? group = ones(Int, n): group
    isempty(y2)? y2 = zeros(Int, n): y2
    isempty(x2)? x2 = zeros(Int, n): x2

    #Concatenate arrays into Dicts
    for i in 1:n
        res[i] = Dict{Any, Any}()
        res[i]["x"] = x[i]
        res[i]["x2"] = x2[i]
        res[i]["y"] = y[i]
        res[i]["group"] = group[i]
        res[i]["y2"] = y2[i]
    end

    #Set VegaData with values, create/add to VegaVisualization
    d = VegaData(name = name, values = res)
    v.data == nothing? v.data = [d]: push!(v.data, d)

    return v
end

max_stars_repo_name: jpthor/vega


text:
#=
    This file includes all the funtions needed to compute the operation
    and investment dynamics
 =#

function compute_operation_dynamics!(h::Int64, y::Int64, s::Int64, mg::Microgrid, controller::AbstractController)
    # Converters
    for (k, a) in enumerate(mg.converters)
        compute_operation_dynamics!(h, y, s, a, controller.decisions.converters[k][h,y,s], mg.parameters.Δh)
    end
    # Storage
    for (k, a) in enumerate(mg.storages)
        compute_operation_dynamics!(h, y, s, a, controller.decisions.storages[k][h,y,s], mg.parameters.Δh)
    end
end
function compute_investment_dynamics!(y::Int64, s::Int64, mg::Microgrid, designer::AbstractDesigner)
     # Generations
     for (k, a) in enumerate(mg.generations)
        compute_investment_dynamics!(y, s, a, designer.decisions.generations[k][y,s])
     end
     # Storages
     for (k, a) in enumerate(mg.storages)
         compute_investment_dynamics!(y, s, a, designer.decisions.storages[k][y,s])
     end
     # Converters
     for (k, a) in enumerate(mg.converters)
         compute_investment_dynamics!(y, s, a, designer.decisions.converters[k][y,s])
     end
end


text:
using Pkg
Pkg.add("Documenter")
using Documenter, QuantizedArrays

# Make src directory available
push!(LOAD_PATH,"../src/")

# Make documentation
makedocs(
    modules = [QuantizedArrays],
    format = Documenter.HTML(),
    sitename = "  ",
    authors = "Corneliu Cofaru, 0x0α Research",
    clean = true,
    debug = true,
    pages = [
        "Introduction" => "index.md",
        "Usage examples" => "examples.md",
        "API Reference" => "api.md",
    ]
)

# Deploy documentation
deploydocs(
    repo = "github.com/zgornel/QuantizedArrays.jl.git",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
module CellArrayApply

using CellwiseValues
using CellwiseValues.CellNumberApply: _checks
using CellwiseValues.CellNumberApply: _getvalues
using CellwiseValues.Kernels: _size_for_broadcast

import CellwiseValues: apply
import Base: iterate
import Base: length
import Base: size
import Base: getindex
import Base: IndexStyle
import CellwiseValues.CellNumberApply: _apply

function apply(k::ArrayKernel,v::Vararg{<:CellValue})
  CellArrayFromKernel(k,v...)
end

function apply(k::ArrayKernel,v::Vararg{<:IndexCellValue})
  IndexCellArrayFromKernel(k,v...)
end

function _apply(f,v,::Val{true})
  k = ArrayKernelFromBroadcastedFunction(f)
  apply(k,v...)
end

struct CellArrayFromKernel{T,N,K,V} <: IterCellArray{T,N,CachedArray{T,N,Array{T,N}}}
  kernel::K
  cellvalues::V
end

function CellArrayFromKernel(k::ArrayKernel,v::Vararg{<:CellValue})
  _checks(v)
  T = _compute_T(k,v)
  N = _compute_N(k,v)
  K = typeof(k)
  V = typeof(v)
  CellArrayFromKernel{T,N,K,V}(k,v)
end

function _compute_T(k,v)
  t = tuple([ _eltype(vi) for vi in v ]...)
  T = compute_type(k,t...)
  @assert T <: NumberLike
  T
end

function _compute_N(k,v)
  d = _compute_ndims(v...)
  compute_ndim(k,d...)
end

_nd(v::CellNumber) = 0

_nd(v::CellArray{T,N}) where {T,N} = N

# TODO use a generated function here
_compute_ndims(v...) = @notimplemented
_compute_ndims(v1) = (_nd(v1),)
_compute_ndims(v1,v2) = (_nd(v1),_nd(v2))
_compute_ndims(v1,v2,v3) = (_nd(v1),_nd(v2),_nd(v3))
_compute_ndims(v1,v2,v3,v4) = (_nd(v1),_nd(v2),_nd(v3),_nd(v4))
_compute_ndims(v1,v2,v3,v4,v5) = (_nd(v1),_nd(v2),_nd(v3),_nd(v4),_nd(v5))
_compute_ndims(v1,v2,v3,v4,v5,v6) = (_nd(v1),_nd(v2),_nd(v3),_nd(v4),_nd(v5),_nd(v6))

_eltype(v::CellNumber{T}) where T = T

_eltype(v::CellArray{T}) where T = T

function length(self::CellArrayFromKernel)
  vi, = self.cellvalues
  length(vi)
end

@inline function iterate(self::CellArrayFromKernel{T,N}) where {T,N}
  zipped = zip(self.cellvalues...)
  znext = iterate(zipped)
  v = CachedArray(T,N)
  _iterate(self,znext,zipped,v)
end

@inline function iterate(self::CellArrayFromKernel,state)
  zipped, zstate, v = state
  znext = iterate(zipped,zstate)
  _iterate(self,znext,zipped,v)
end

@inline function _iterate(self::CellArrayFromKernel,znext,zipped,v)
  if znext === nothing; return nothing end
  a, zstate = znext
  s = _compute_sizes(a...)
  z = compute_size(self.kernel,s...)
  setsize!(v,z)
  compute_value!(v,self.kernel,a...)
  state = (zipped,zstate,v)
  (v, state)
end

const _bs = _size_for_broadcast

# TODO use a generated function?
_compute_sizes(a...) = @notimplemented
_compute_sizes(a1) = (_bs(a1),)
_compute_sizes(a1,a2) = (_bs(a1),_bs(a2))
_compute_sizes(a1,a2,a3) = (_bs(a1),_bs(a2),_bs(a3))
_compute_sizes(a1,a2,a3,a4) = (_bs(a1),_bs(a2),_bs(a3),_bs(a4))
_compute_sizes(a1,a2,a3,a4,a5) = (_bs(a1),_bs(a2),_bs(a3),_bs(a4),_bs(a5))
_compute_sizes(a1,a2,a3,a4,a5,a6) = (_bs(a1),_bs(a2),_bs(a3),_bs(a4),_bs(a5),_bs(a6))

struct IndexCellArrayFromKernel{T,N,K,V} <: IndexCellArray{T,N,CachedArray{T,N,Array{T,N}},1}
  kernel::K
  cellvalues::V
  cache::CachedArray{T,N,Array{T,N}}
end

function IndexCellArrayFromKernel(k::ArrayKernel,v::Vararg{<:CellValue})
  _checks(v)
  T = _compute_T(k,v)
  N = _compute_N(k,v)
  K = typeof(k)
  V = typeof(v)
  cache = CachedArray(T,N)
  IndexCellArrayFromKernel{T,N,K,V}(k,v,cache)
end

function length(self::IndexCellArrayFromKernel)
  vi, = self.cellvalues
  length(vi)
end

size(self::IndexCellArrayFromKernel) = (length(self),)

function getindex(self::IndexCellArrayFromKernel,i::Integer)
  a = _getvalues(i,self.cellvalues...)
  s = _compute_sizes(a...)
  z = compute_size(self.kernel,s...)
  v = self.cache
  setsize!(v,z)
  compute_value!(v,self.kernel,a...)
  v
end

end # module CellArrayApply


text:
# This file was generated, do not modify it. # hide
fitted_params(rfr_tm).best_model


text:
"""
Kinetic energy: 1/2 sum_n f_n ∫ |∇ψn|^2.
"""
struct Kinetic
    scaling_factor::Real
end
Kinetic(; scaling_factor=1) = Kinetic(scaling_factor)
(kin::Kinetic)(basis) = TermKinetic(basis, kin.scaling_factor)
function Base.show(io::IO, kin::Kinetic)
    fac = isone(kin.scaling_factor) ? "" : ", scaling_factor=$scaling_factor"
    print(io, "Kinetic($fac)")
end

struct TermKinetic <: Term
    scaling_factor::Real  # scaling factor, absorbed into kinetic_energies
    kinetic_energies::Vector{<:AbstractVector}  # kinetic energy 1/2|G+k|^2 for every kpoint
end
function TermKinetic(basis::PlaneWaveBasis{T}, scaling_factor) where {T}
    kinetic_energies = [[T(scaling_factor) * sum(abs2, Gk) / 2
                         for Gk in Gplusk_vectors_cart(basis, kpt)]
                        for kpt in basis.kpoints]
    TermKinetic(T(scaling_factor), kinetic_energies)
end

@timing "ene_ops: kinetic" function ene_ops(term::TermKinetic, basis::PlaneWaveBasis{T},
                                            ψ, occ; kwargs...) where {T}
    ops = [FourierMultiplication(basis, kpoint, term.kinetic_energies[ik])
           for (ik, kpoint) in enumerate(basis.kpoints)]
    isnothing(ψ) && return (E=T(Inf), ops=ops)

    E = zero(T)
    for (ik, k) in enumerate(basis.kpoints)
        for iband = 1:size(ψ[1], 2)
            ψnk = @views ψ[ik][:, iband]
            E += (basis.kweights[ik] * occ[ik][iband]
                  * real(dot(ψnk, Diagonal(term.kinetic_energies[ik]), ψnk)))
        end
    end
    E = mpi_sum(E, basis.comm_kpts)

    (E=E, ops=ops)
end


text:
function generate_eval_expr(domain::Domain, state::State, term::Term,
                            varmap=Dict{Var,Any}(), state_var=:state)
    if (term.name in keys(get_fluents(domain)) || term isa Var)
        return generate_get_expr(domain, state, term, varmap, state_var)
    elseif term isa Const
        return QuoteNode(term.name)
    end
    subexprs = [generate_eval_expr(domain, state, a, varmap, state_var)
                for a in term.args]
    expr = if is_global_func(term)
        op = GLOBAL_FUNCTIONS[term.name]
        Expr(:call, QuoteNode(op), subexprs...)
    else
        error("Unrecognized predicate or operator $(term.name).")
    end
    return expr
end

function generate_quantified_expr(domain::Domain, state::State, term::Term,
                                  varmap=Dict{Var,Any}(), state_var=:state)
    @assert length(term.args) == 2 "$(term.name) takes two arguments"
    varmap = copy(varmap) # Make local copy of variable context
    typeconds, query = flatten_conjs(term.args[1]), term.args[2]
    types, vars = Symbol[], Symbol[]
    for (i, cond) in enumerate(typeconds)
        v = Symbol("o$i") # Create local variable name
        push!(vars, v)
        push!(types, cond.name) # Extract type name
        varmap[cond.args[1]] = :($v.name)
    end
    # Generate query expression with local variable context
    expr = generate_check_expr(domain, state, query, varmap, state_var)
    # Construct iterator over (typed) objects
    v, ty = pop!(vars), QuoteNode(pop!(types))
    expr = :($expr for $v in get_objects(state, $ty))
    while !isempty(types)
        v, ty = pop!(vars), QuoteNode(pop!(types))
        expr = Expr(:flatten, :($expr for $v in get_objects($state_var, $ty)))
    end
    if Domain isa AbstractedDomain
        accum = term.name == :forall ? :(&) : :(|)
        return :($accum($expr...))
    else
        accum = term.name == :forall ? :all : :any
        return :($accum($expr))
    end
end

function generate_check_expr(domain::Domain, state::State, term::Term,
                             varmap=Dict{Var,Any}(), state_var=:state)
    if (is_global_func(term) || term.name in keys(get_funcdefs(domain)))
        return generate_eval_expr(domain, state, term, varmap, state_var)
    elseif (term.name in keys(get_fluents(domain)) || term isa Var)
        return generate_get_expr(domain, state, term, varmap, state_var)
    elseif term.name in (:forall, :exists)
        return generate_quantified_expr(domain, state, term, varmap, state_var)
    elseif term isa Const
        return QuoteNode(term.name::Bool)
    end
    subexprs = [generate_check_expr(domain, state, a, varmap, state_var)
                for a in term.args]
    expr = if term.name == :and
        foldr((a, b) -> Expr(:&&, a, b), subexprs)
    elseif term.name == :or
        foldr((a, b) -> Expr(:||, a, b), subexprs)
    elseif term.name == :imply
        @assert length(term.args) == 2 "imply takes two arguments"
        :(!$(subexprs[1]) || $(subexprs[2]))
    elseif term.name == :not
        @assert length(term.args) == 1 "not takes one argument"
        :(!$(subexprs[1]))
    else
        error("Unrecognized predicate or operator $(term.name).")
    end
    return expr
end

function generate_check_expr(domain::AbstractedDomain, state::State, term::Term,
                             varmap=Dict{Var,Any}(), state_var=:state)
    if (is_global_func(term) || term.name in keys(get_funcdefs(domain)))
        return generate_eval_expr(domain, state, term, varmap, state_var)
    elseif (term.name in keys(get_fluents(domain)) || term isa Var)
        return generate_get_expr(domain, state, term, varmap, state_var)
    elseif term.name in (:forall, :exists)
        return generate_quantified_expr(domain, state, term, varmap, state_var)
    elseif term isa Const
        return QuoteNode(term.name::Bool)
    end
    subexprs = [generate_check_expr(domain, state, a, varmap, state_var)
                for a in term.args]
    expr = if term.name == :and
        foldr((a, b) -> Expr(:call, :&, a, b), subexprs)
    elseif term.name == :or
        foldr((a, b) -> Expr(:call, :|, a, b), subexprs)
    elseif term.name == :imply
        @assert length(term.args) == 2 "imply takes two arguments"
        :(!$(subexprs[1]) | $(subexprs[2]))
    elseif term.name == :not
        @assert length(term.args) == 1 "not takes one argument"
        :(!$(subexprs[1]))
    else
        error("Unrecognized predicate or operator $(term.name).")
    end
    return expr
end

function generate_forall_effect_expr(domain::Domain, state::State, term::Term,
                                     varmap=Dict{Var,Any}(), state_var=:state)
    @assert length(term.args) == 2 "$(term.name) takes two arguments"
    varmap = copy(varmap) # Make local copy of variable context
    typeconds, effect = flatten_conjs(term.args[1]), term.args[2]
    types, vars = Symbol[], Symbol[]
    for (i, cond) in enumerate(typeconds)
        v = Symbol("o$i") # Create local variable name
        push!(vars, v)
        push!(types, cond.name) # Extract type name
        varmap[cond.args[1]] = :($v.name)
    end
    # Generate effect expression with local variable context
    expr = generate_effect_expr(domain, state, effect, varmap, state_var)
    # Special case if only one object variable is enumerated over
    if length(vars) == 1
        v, ty = vars[1], QuoteNode(types[1])
        return quote for $v in get_objects($state_var, $ty)
            $expr
        end end
    end
    # Construct iterator over (typed) objects
    obj_exprs = (:(get_objects($state_var, $(QuoteNode(ty)))) for ty in types)
    expr = quote for ($(vars...),) in zip($(obj_exprs...))
        $expr
    end end
    return expr
end

function generate_effect_expr(domain::Domain, state::State, term::Term,
                              varmap=Dict{Var,Any}(), state_var=:state)
    expr = if term.name == :and
        subexprs = [generate_effect_expr(domain, state, a, varmap, state_var)
                    for a in term.args]
        Expr(:block, subexprs...)
    elseif term.name == :assign
        @assert length(term.args) == 2 "assign takes two arguments"
        term, val = term.args
        val = generate_eval_expr(domain, state, val, varmap, :prev_state)
        generate_set_expr(domain, state, term, val, varmap, state_var)
    elseif term.name == :not
        @assert length(term.args) == 1 "not takes one argument"
        term = term.args[1]
        @assert term.name in keys(get_predicates(domain)) "unrecognized predicate"
        generate_set_expr(domain, state, term, false, varmap, state_var)
    elseif term.name == :when
        @assert length(term.args) == 2 "when takes two arguments"
        cond, effect = term.args
        cond = generate_check_expr(domain, state, cond, varmap, :prev_state)
        effect = generate_effect_expr(domain, state, effect, varmap, state_var)
        :(cond = $cond; (cond == true || cond == both) && $effect)
    elseif term.name == :forall
        generate_forall_effect_expr(domain, state, term, varmap, state_var)
    elseif term.name in keys(GLOBAL_MODIFIERS)
        @assert length(term.args) == 2 "$(term.name) takes two arguments"
        op = GLOBAL_MODIFIERS[term.name]
        term, val = term.args
        prev_val = generate_get_expr(domain, state, term, varmap, :prev_state)
        val = generate_eval_expr(domain, state, val, varmap, :prev_state)
        new_val = Expr(:call, op, prev_val, val)
        generate_set_expr(domain, state, term, new_val, varmap, state_var)
    elseif term.name in keys(get_predicates(domain))
        generate_set_expr(domain, state, term, true, varmap, state_var)
    else
        error("Unrecognized predicate or operator $(term.name).")
    end
    return expr
end


text:
#=
epigraph of perspective of a separable spectral function
on a real vector or real symmetric or complex Hermitian (svec scaled triangle) domain
(u, v, w) : v ≥ 0, W ⪰ 0, u ≥ v * h(W / v)

dual cone is
(u, v, w) : u ≥ 0, W ∈ dom(h⋆), v ≥ u * h⋆(W / u)
=#

import Hypatia.Cones: VectorCSqr, MatrixCSqr, SepSpectralFun
import Hypatia.Cones: h_val, h_conj_dom_pos, h_conj, h_der1

include("vector.jl")
include("matrix.jl")

function MOIPajarito.Cones.create_cache(
    oa_s::Vector{AE},
    cone::Hypatia.EpiPerSepSpectralCone{RealF},
    opt::Optimizer,
)
    @assert MOI.dimension(cone) == length(oa_s)
    cache = create_sepspectral_cache(cone.Q, cone.use_dual, cone.d, opt)
    cache.oa_s = oa_s
    cache.d = cone.d
    cache.h = cone.h
    return cache
end

function per_sepspec(f::Function, h::SepSpectralFun, v::RealF, w::AbstractVector{RealF})
    v < 1e-12 && return 0.0 # TODO?
    if v == 1
        return f(w, h)
    end
    return v * f(w / v, h)
end

swap_epiper(::Type{Prim}, u::T, v::T) where {T} = (u, v)
swap_epiper(::Type{Dual}, u::T, v::T) where {T} = (v, u)

dom_pos(::Type{Prim}, ::SepSpectralFun) = true
dom_pos(::Type{Dual}, h::SepSpectralFun) = h_conj_dom_pos(h)

val_or_conj(::Type{Prim}) = h_val
val_or_conj(::Type{Dual}) = h_conj
conj_or_val(::Type{Prim}) = h_conj
conj_or_val(::Type{Dual}) = h_val

import Hypatia.Cones: NegLogSSF, NegEntropySSF, NegSqrtSSF, NegPower01SSF, Power12SSF

# r values for initial cuts
init_r_vals(::Type{<:Union{Prim, Dual}}, ::NegLogSSF) = [1e-4, 1.0, 1e4]
init_r_vals(::Type{Prim}, ::NegEntropySSF) = [-9.0, -1.0, 8.0]
init_r_vals(::Type{Dual}, ::NegEntropySSF) = [1e-4, 1.0, 1e3]
init_r_vals(::Type{<:Union{Prim, Dual}}, ::NegSqrtSSF) = [1e-4, 1.0, 1e4]

function init_r_vals(::Type{Prim}, h::NegPower01SSF)
    vals = [1e-4, 1.0, 1e4]
    p = RealF(h.p)
    p <= 0.5 && return vals
    q = p / (p - 1)
    return [p * (v / (1 - p))^inv(q) for v in vals]
end

init_r_vals(::Type{Dual}, h::NegPower01SSF) = [1e-4, 1.0, 1e4]

function init_r_vals(::Type{Prim}, h::Power12SSF)
    p = RealF(h.p)
    p >= 1.5 && return [-1e2, -3.0, 1.0]
    q = p / (p - 1)
    vals = [-1e4, -1.0, 1.0]
    return [sign(v) * p * (abs(v) / (p - 1))^inv(q) for v in vals]
end

init_r_vals(::Type{Dual}, ::Power12SSF) = [1e-3, 1.0, 1e3]

# gradients
h_grad(::Type{Prim}, h::SepSpectralFun, w::Vector{RealF}) = h_der1(similar(w), w, h)

h_grad(::Type{Dual}, h::NegLogSSF, w::Vector{RealF}) = h_der1(similar(w), w, h)
h_grad(::Type{Dual}, ::NegEntropySSF, w::Vector{RealF}) = [-exp(-w_i - 1) for w_i in w]
h_grad(::Type{Dual}, ::NegSqrtSSF, w::Vector{RealF}) = [-0.25 * w_i^-2 for w_i in w]

function h_grad(::Type{Dual}, h::NegPower01SSF, w::Vector{RealF})
    p = RealF(h.p)
    qm1 = inv(p - 1)
    c = -p^-qm1
    return [c * w_i^qm1 for w_i in w]
end

function h_grad(::Type{Dual}, h::Power12SSF, w::Vector{RealF})
    p = RealF(h.p)
    qm1 = inv(p - 1)
    c = -p^-qm1
    g = zero(w)
    for (i, w_i) in enumerate(w)
        if w_i < 0
            g[i] = c * abs(w_i)^qm1
        end
    end
    return g
end


text:
using SmartShockFinder

datatyp = "NNRH"
data_input_size = 5
n_meshes = 4
polydeg = collect(2:9)

generate_validdataset1d(datatyp, n_meshes, polydeg, data_input_size)


text:
# functions to determine hermite coefficients of a function

hermite(::Type{Val{0}}) = 1.0
hermite(::Type{Val{0}}, ξ) = hermite(Val{0})
hermite(::Type{Val{1}}, ξ) = ξ .* hermite(Val{0})

hermite(::Type{Val{0}}, ξ, q) = hermite(Val{0})
hermite(::Type{Val{1}}, ξ, q) = ξ .* hermite(Val{0})

hermite(N::Int, ξ) = hermite(Val{N}, ξ)
hermite(N::Int, ξ, q::Quadrature) = hermite(Val{N}, ξ, q)

function hermite(::Type{Val{2}}, ξ)
    D = length(ξ)
    return [ξ[β] * ξ[α] - δ(α, β) for α in 1:D, β in 1:D]
end
function hermite(::Type{Val{3}}, ξ)
    D = length(ξ)
    return [
        ξ[γ] * ξ[β] * ξ[α] - (ξ[α] * δ(β, γ) + ξ[β] * δ(α, γ) + ξ[γ] * δ(α, β))

        for α in 1:D, β in 1:D, γ in 1:D
    ]
end

function hermite(::Type{Val{4}}, ξ)
    D = length(ξ)
    return [
        ξ[ϵ] * ξ[γ] * ξ[β] * ξ[α] - (ξ[α] *
        ξ[β] *
        δ(γ, ϵ) + ξ[α] *
        ξ[γ] *
        δ(β, ϵ) + ξ[α] *
        ξ[ϵ] *
        δ(β, γ) + ξ[β] *
        ξ[γ] *
        δ(α, ϵ) + ξ[β] *
        ξ[ϵ] *
        δ(α, γ) + ξ[γ] *
        ξ[ϵ] *
        δ(α, β)) + (δ(α, β) * δ(γ, ϵ) + δ(α, γ) * δ(β, ϵ) + δ(α, ϵ) * δ(β, γ))

        for α in 1:D, β in 1:D, γ in 1:D, ϵ in 1:D
    ]
end

function hermite(::Type{Val{2}}, ξ, q)
    D = length(ξ)
    cs = 1 / q.speed_of_sound_squared
    return [ξ[β] * ξ[α] - cs * δ(α, β) for α in 1:D, β in 1:D]
end
function hermite(::Type{Val{3}}, ξ, q)
    D = length(ξ)
    cs = 1 / q.speed_of_sound_squared
    return [
        ξ[γ] * ξ[β] * ξ[α] - cs * (ξ[α] * δ(β, γ) + ξ[β] * δ(α, γ) + ξ[γ] * δ(α, β))

        for α in 1:D, β in 1:D, γ in 1:D
    ]
end

function hermite(::Type{Val{4}}, ξ, q)
    D = length(ξ)
    cs = 1 / q.speed_of_sound_squared
    return [
        ξ[ϵ] * ξ[γ] * ξ[β] * ξ[α] - cs * (ξ[α] *
        ξ[β] *
        δ(γ, ϵ) + ξ[α] *
        ξ[γ] *
        δ(β, ϵ) + ξ[α] *
        ξ[ϵ] *
        δ(β, γ) + ξ[β] *
        ξ[γ] *
        δ(α, ϵ) + ξ[β] *
        ξ[ϵ] *
        δ(α, γ) + ξ[γ] *
        ξ[ϵ] *
        δ(α, β)) + cs^2 * (δ(α, β) * δ(γ, ϵ) + δ(α, γ) * δ(β, ϵ) + δ(α, ϵ) * δ(β, γ))

        for α in 1:D, β in 1:D, γ in 1:D, ϵ in 1:D
    ]
end


text:
type SGksd <: StepData
	getStep::Function
	stepString::AbstractString
	function SGksd(stepString::AbstractString)
		new(computeSGkStep,stepString)
	end
	function SGksd()
		new(computeSGkStep,"SGk")
	end
end

function computeSGkStep(x, k, gnum, sd::SGksd, problem::Problem; outputLevel =0)

	func = consume(problem.getNextSampleFunction)

	(f, g) = func(x)

	gnum += 1

	(g/(k+1), gnum)

end


text:
#supplybounce.jl: Generate plots of supply bounce
#-------------------------------------------------------------------------------

c2rgba(c::Colorant) = convert(Colors.RGBA{Colors.FixedPointNumbers.N0f8}, c)

#Perform color mix:
function cmix(c1::T, c2::T, w1::Float64) where T <: Colors.RGBA{Colors.FixedPointNumbers.N0f8}
	w2 = 1-w1
	return T(c1.r*w1+c2.r*w2, c1.g*w1+c2.g*w2, c1.b*w1+c2.b*w2, c1.alpha*w1+c2.alpha*w2)
end

cmix(c1::Colorant, c2::Colorant, w1::Float64) =
	cmix(c2rgba(c1), c2rgba(c2), w1)

#==Constants
===============================================================================#
line_sup = cons(:a, line = set(style=:solid, color=:red, width=3))
line_gnd = cons(:a, line = set(style=:solid, color=:blue, width=3))
line_sig = cons(:a, line = set(style=:solid, color=:green, width=3))

line_idealsup = cons(:a, line = set(style=:dash, color=cmix(colorant"red", colorant"white", .5), width=2))
line_idealgnd = cons(:a, line = set(style=:dash, color=cmix(colorant"blue", colorant"white", .5), width=2))

#==
===============================================================================#
function plot_supplybounce(plotname::String="supplybounce";
		datapath=joinpath(@__DIR__, "..", "..", "data", "SupplyTransients.csv")
	)
	data = readdlm(datapath, ',', comments=true)
	#Time, VDD, VSS, VSIG

	t = data[:, 1]
	VDD = DataF1(t, data[:, 2])
	VSS = DataF1(t, data[:, 3])
	Vload = DataF1(t, data[:, 4])
	VDD_diff = VDD-VSS
	simgnd = 0*VDD
	VDD_ideal = simgnd+1

	#NOTE: Using 1.25 to avoid bug displaying 1000e-3V
	plotattr = cons(:a,
		ystrip1 = set(scale=:lin, axislabel=LBLAX_VOLT, min=-0.25, max=1.25),
		ystrip2 = set(scale=:lin, axislabel=LBLAX_VOLT, min=-0.25, max=1.25),
		ystrip3 = set(scale=:lin, axislabel=LBLAX_VOLT, min=-0.25, max=1.25),
#		xaxis = set(max = 10e-9),
	)

	tmid = 10e-9
	midpos = (tmid,0.5)
	siglist=[
		(id="(single Lpkg)", wfrm=VDD_ideal, line=line_idealsup, idpos=midpos, posalign=:cc, strip=2),
		(id="(LVDD=LVSS=Lpkg/2)", wfrm=VDD_ideal, line=line_idealsup, idpos=midpos, posalign=:cc, strip=3),
		(id="", wfrm=simgnd, line=line_idealgnd, idpos=(0,0), posalign=:bl, strip=2),
		(id="", wfrm=simgnd, line=line_idealgnd, idpos=(0,0), posalign=:bl, strip=3),

		(id="Vload", wfrm=Vload, line=line_sig, idpos=midpos, posalign=:cl, strip=1),
		(id="VDD", wfrm=VDD_diff, line=line_sup, idpos=(tmid, 0.9), posalign=:tc, strip=2),
		(id="VDD", wfrm=VDD, line=line_sup, idpos=(tmid,0.95), posalign=:tc, strip=3),
		(id="VSS", wfrm=VSS, line=line_gnd, idpos=(tmid,0.05), posalign=:bc, strip=3),
	]

	data = (
		nstrips = 3,
		plotattr,
		siglist=siglist,
	)

	plot = EasyPlot.build(plot_builders[:tran], data)
	plotgui = display(plotdisplay[:dflt], plot)
	EasyPlot._write(:png, "$plotname.png", optnoleg, plotgui)
	EasyPlot._write(:eps, "$plotname.eps", optnoleg, plotgui)
	return plotgui
end

#last line


text:
#############################################################################
# JuMP
# An algebraic modeling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# mindistortion.jl
#
# This example arises from computational geometry, in particular the problem
# of embedding a general finite metric space into a euclidean space.
#
# It is known that the 4-point metric space defined by the star graph:
#   x
#    \
#     x — x
#    /
#   x
# where distances are computed by length of the shortest path between
# vertices, cannot be exactly embedded into a euclidean space of any
# dimension.
#
# Here we will formulate and solve an SDP to compute the best possible
# embedding, that is, the embedding f() that minimizes the distortion c such
# that (1/c)*D(a,b) ≤ ||f(a)-f(b)|| ≤ D(a,b) for all points (a,b), where
# D(a,b) is the distance in the metric space.
#
# Any embedding can be characterized by its Gram matrix Q, which is PSD,
# and ||f(a)-f(b)||^2 = Q[a,a] + Q[b,b] - 2Q[a,b]
# We can therefore constrain
# D[i,j]^2 ≤ Q[i,i] + Q[j,j] - 2Q[i,j] ≤ c^2*D[i,j]^2
# and minimize c^2, which gives us the SDP formulation below.
#
# For more detail, see
# "Lectures on discrete geometry" by J. Matoušek, Springer, 2002, pp. 378-379.

using JuMP

m = Model()

D = [0.0 1.0 1.0 1.0
     1.0 0.0 2.0 2.0
     1.0 2.0 0.0 2.0
     1.0 2.0 2.0 0.0]

@variable(m, cSq >= 1.0)

@variable(m, Q[1:4,1:4], SDP)

for i in 1:4
    for j in (i+1):4
        @constraint(m, D[i,j]^2 <= Q[i,i] + Q[j,j] - 2Q[i,j])
        @constraint(m, Q[i,i] + Q[j,j] - 2Q[i,j] <= D[i,j]^2*cSq )
    end
end

@objective(m, Min, cSq)

solve(m)

println(getvalue(cSq))


text:
include(find_in_path("PatternDispatch.jl"))

module TestShowPattern
using PatternDispatch

@show (@qpat 1)
@show (@qpat 1::Int)
@show (@qpat 1::Float64)

println()
@show (@qpat x::Any)
@show (@qpat ::Any)
@show (@qpat ::None)
@show (@qpat (x::None,y::Int))

println()
@show (@qpat x)
@show (@qpat ::Int)
@show (@qpat x::Int)
@show (@qpat x~y)
@show (@qpat x~y::Int)

println()
@show (@qpat (x,y))
@show (@qpat (x,y)::Tuple)
@show (@qpat z~(x::Int,1))

end


text:
module Apriori

using DataFrames
using Combinatorics
using CSV

ATTRIBUTE_TYPE = Int32

function apriori(
    data::DataFrame,
    min_relative_support=0.2,
    min_confidence=0.3;
    X_in_antecedent=Vector{Symbol}()::Vector{Symbol},
    antecedent_in_X=Vector{Symbol}()::Vector{Symbol},
    X_in_consequent=Vector{Symbol}()::Vector{Symbol},
    consequent_in_X=Vector{Symbol}()::Vector{Symbol},
    )::Vector{Pair{Tuple{Vector{Symbol},Vector{Symbol}},Tuple{Float64,Float64}}}

    attributes = propertynames(data)::Array{Symbol,1}

    symbol2index = Dict([x => y for (x, y) in zip(attributes, one(ATTRIBUTE_TYPE):convert(ATTRIBUTE_TYPE, length(attributes)))])

    filters = __prepare_filters(symbol2index, X_in_antecedent, antecedent_in_X, X_in_consequent, consequent_in_X)

    frequent_itemsets = __apriori_frequent_itemsets(data, min_relative_support)

    rules = __apriori_rule_gen(frequent_itemsets, min_confidence, filters)
    map(x -> translate_rule(x, attributes), rules)
end

function apriori_frequent_itemsets(data::DataFrame, min_relative_support=0.2)
    attributes = propertynames(data)
    freq_itemsets = __apriori_frequent_itemsets(data, min_relative_support)
    symbol2index = Dict([x => y for (x, y) in zip(attributes, 1:ncol(data))])
    frequent_itemsets = map(x -> Set{Symbol}(getindex.(Ref(attributes), x[1])) => x[2], freq_itemsets)
end

function __apriori_frequent_itemsets(data::DataFrame, min_relative_support=0.2)::Vector{Pair{Vector{ATTRIBUTE_TYPE},Int64}}
    supp = floor(min_relative_support * nrow(data))

    frequent_itemsets = Vector{Pair{Vector{ATTRIBUTE_TYPE},Int64}}([[]=>nrow(data)])

    itemsets = [Vector{ATTRIBUTE_TYPE}([i]) for i in 1:ncol(data)]
    infrequents = Set{Vector{ATTRIBUTE_TYPE}}()

    while true
        itesmets_w_support = map(x -> x => support(data, x), itemsets)
        frequent_itesmets_w_support = filter(x -> x[2] > supp, itesmets_w_support)
        append!(frequent_itemsets, frequent_itesmets_w_support)

        if (length(frequent_itesmets_w_support) <= 1) break end

        itemsets = filter(x -> !any(subsets(x, length(x) - 1) .∈ (infrequents,)), merge_vectors(map(x -> x[1], frequent_itesmets_w_support)))
        infrequents = map(y -> y[1], filter(x -> x[2] <= supp, itesmets_w_support))
    end

    frequent_itemsets
end

function __prepare_filters(translate_dict,
    X_in_antecedent=Vector{Symbol}()::Vector{Symbol},
    antecedent_in_X=Vector{Symbol}()::Vector{Symbol},
    X_in_consequent=Vector{Symbol}()::Vector{Symbol},
    consequent_in_X=Vector{Symbol}()::Vector{Symbol},
    )

    t = [
        :x_in_ant => X_in_antecedent,
        :ant_in_x => antecedent_in_X,
        :x_in_con => X_in_consequent,
        :con_in_x => consequent_in_X
    ]

    @assert all(X_in_antecedent .∉ (X_in_consequent,))

    filters = Dict(filter(x -> length(x[2]) > 0, map(x -> x[1] => sort(translate_itemset(x[2], translate_dict)), t)))
end

function apriori_rule_gen(
    frequent_itemsets::Vector{Pair{Set{Symbol},supp_int}},
    attribute_names::Vector{Symbol},
    min_confidence=0.3;
    X_in_antecedent=Vector{Symbol}()::Vector{Symbol},
    antecedent_in_X=Vector{Symbol}()::Vector{Symbol},
    X_in_consequent=Vector{Symbol}()::Vector{Symbol},
    consequent_in_X=Vector{Symbol}()::Vector{Symbol},
    ) where {supp_int <: Signed}

    symbol2index = Dict([x => convert(ATTRIBUTE_TYPE, y) for (x, y) in zip(attribute_names, 1:length(attribute_names))])
    frequent_itemsets = map(x -> sort(Vector{ATTRIBUTE_TYPE}(getindex.(Ref(symbol2index), x[1]))) => x[2], frequent_itemsets)

    filters = __prepare_filters(symbol2index, X_in_antecedent, antecedent_in_X, X_in_consequent, consequent_in_X)

    rules = __apriori_rule_gen(frequent_itemsets, min_confidence, filters)
    result = map(x -> translate_rule(x, attribute_names), rules)
    result::Vector{Pair{Tuple{Vector{Symbol},Vector{Symbol}},Float64}}
end

function __apriori_rule_gen(
    frequent_itemsets::Vector{Pair{Vector{ATTRIBUTE_TYPE},supp_int}},
    min_confidence=0.3,
    filters=Dict()
    ) where {supp_int <: Signed}

    df = Dict(frequent_itemsets)

    dataset_size = df[[]]

    all_strong_rules = Vector{Pair{Tuple{Vector{ATTRIBUTE_TYPE},Vector{ATTRIBUTE_TYPE}},Tuple{Float64,Float64}}}()

    for (Z_set, Z_sup) in frequent_itemsets
        #if (length(Z_set) == 1) continue end

        Z = Z_set
        i = 1
        Y = subsets(Z, i)

        if (haskey(filters, :x_in_ant))
            X = filters[:x_in_ant]
            if (all(X .∈ (Z,)))
                Y = filter(x -> x[end] ∉ X, Y)
            else continue end
        end

        if (haskey(filters, :con_in_x))
            X = filters[:con_in_x]
            Y = filter(x -> x[end] ∈ X, Y)
        end

        if (haskey(filters, :ant_in_x))
            X = filters[:ant_in_x]
            X = filter(x -> x[1] ∈ Z, X)
            ZmX = filter(x -> x ∉ X, Z)
            if (length(ZmX) > 0 && length(X) > 0)
                sp = ((X, ZmX) => (Z_sup/dataset_size, Z_sup / df[X]))
                if (sp[2][2] >= min_confidence &&
                    (!haskey(filters, :x_in_ant) || all(filters[:x_in_ant] .∈ (sp[1][1],))) && # ghetto fix na inne filtry żeby obejmowały te wyjątki, ale nadal są powtórzenia wśród tych special cases
                    (!haskey(filters, :con_in_x) || all(sp[1][2] .∈ (filters[:con_in_x],))) &&
                    (!haskey(filters, :x_in_con) || all(filters[:x_in_con] .∈ (sp[1][2],))) &&
                    sp ∉ all_strong_rules)
                    append!(all_strong_rules,  [sp])
                end
            end
            Y = filter(x -> x[end] ∈ X, Y)
            Y = map(x -> vcat(ZmX, x), Y)
        end

        if (haskey(filters, :x_in_con))
            X = filters[:x_in_con]
            if (all(X .∈ (Z,))) # czy wszystkie elementy X są w Z
                ZmX = filter(x -> x ∉ X, Z)
                if (length(ZmX) > 0)
                    sp = (ZmX, X) => (Z_sup/dataset_size, Z_sup / df[ZmX] )
                    if (sp[2][2] >= min_confidence &&
                        (!haskey(filters, :x_in_ant) || all(filters[:x_in_ant] .∈ (sp[1][1],))) && # ghetto solution
                        (!haskey(filters, :con_in_x) || all(sp[1][2] .∈ (filters[:con_in_x],))) &&
                        (!haskey(filters, :ant_in_x) || all(sp[1][1] .∈ (filters[:ant_in_x],))) &&
                        sp ∉ all_strong_rules)
                        append!(all_strong_rules, [sp])
                    end
                end
                Y = filter(x -> x[end] ∉ X, Y)
                Y = map(x -> vcat(X, x), Y)
            else continue end
        end

        if (length(Y) > 0 &&  all(length.(Y) .< length(Z)))
            i = length(Y[1])
            while true
                strong_rules = filter(x -> x[2][2] >= min_confidence, map(x -> begin
                    ant = antecedent(Z, x)
                    (ant, x) => (Z_sup/dataset_size, Z_sup / df[ant] )
                end, Y))

                append!(all_strong_rules, strong_rules)
                i += 1
                if (i == length(Z_set) || length(strong_rules) == 0) break end

                Y = merge_vectors(map(x -> x[1][2], strong_rules))
            end
        end
    end
    all_strong_rules::Vector{Pair{Tuple{Vector{ATTRIBUTE_TYPE},Vector{ATTRIBUTE_TYPE}},Tuple{Float64,Float64}}}
end

function merge_vectors(itemsets::Vector{Vector{ATTRIBUTE_TYPE}})
    result = Vector{Vector{ATTRIBUTE_TYPE}}()
    itemsets_len = length(itemsets)::Int64
    for i in 1:itemsets_len
        itemset1 = itemsets[i]
        prefix = itemset1[1:end - 1]
        suffixes = Vector{Vector{ATTRIBUTE_TYPE}}()
        for j in i + 1:itemsets_len
            itemset2 = itemsets[j]
            if prefix != itemset2[1:end - 1]
                break
            end
            push!(suffixes, sort([itemset1[end], itemset2[end]]))
        end
        merged_itemsets = map(x -> vcat(prefix, x), sort(suffixes))
        append!(result, merged_itemsets)
    end
    result::Vector{Vector{ATTRIBUTE_TYPE}}
end

function translate_itemset(element::Vector{Symbol}, dict::Dict{Symbol,ATTRIBUTE_TYPE})
    Vector{ATTRIBUTE_TYPE}(getindex.(Ref(dict), element))
end

function translate_rule(rule::Pair{Tuple{Vector{ATTRIBUTE_TYPE},Vector{ATTRIBUTE_TYPE}},Tuple{Float64,Float64}}, dict::Vector{Symbol})
    (getindex.(Ref(dict), rule[1][1]), getindex.(Ref(dict), rule[1][2])) => rule[2]
end

function support(data::DataFrame, attrs::Vector{ATTRIBUTE_TYPE})
    nrow(data[reduce((x, y) -> x .& y, [data[!, attr] .== true for attr in attrs]),:])
end

function subsets(vec::Vector{ATTRIBUTE_TYPE}, len)
    Vector.(combinations([vec...], len))
end

function antecedent(fvec, vec)
    filter(x -> x ∉ vec, fvec)
end

function dummy_dataset(attrs, rows)
    @assert attrs <= 26 "TODO the attr label fun to generate more"
    attr_names = collect('a':'z')[1:attrs]
    DataFrame([Symbol(i) => rand(Bool, rows) for i in attr_names])
end

# Creates a semi-random dataset that exhibits a strong association rule
function dummy_dataset_biased(attrs, rows, relative_sup, antecedent, consequent, support, confidence=1)
    @assert attrs <= 26 "TODO the attr label fun to generate more"
    attr_names = collect('a':'z')[1:attrs]
    df = DataFrame([Symbol(i) => map(x->x<relative_sup, rand(rows)) for i in attr_names])
    # print(df)
    antecedent_indices = findall(col_name -> col_name in antecedent, attr_names)
    consequent_indices = findall(col_name -> col_name in consequent, attr_names)
    antecedent_count = nrow(filter(row -> all([row[index] for index in antecedent_indices]), df))
    rule_count = nrow(filter(row -> all([row[index] for index in vcat(antecedent_indices, consequent_indices)]), df))
    required_rule_count = rows * support

    remove_consequent_count = required_rule_count * (1-confidence)
    for row in eachrow(df)
        if !all([row[index] for index in vcat(antecedent_indices, consequent_indices)])
            for index in vcat(antecedent_indices, consequent_indices)
                row[index] = true
            end
            if remove_consequent_count>0
                row[consequent_indices[1]] = false
                remove_consequent_count -= 1
            end
            rule_count = rule_count + 1
        end
        if rule_count >= required_rule_count
            break
        end
    end
    return df
end

function mushroom_dataset()
    download("https://gist.githubusercontent.com/krynju/120e623b7b0277e45126baf8aa39edad/raw/f980a374b787cd7e539e821450e4de91f948b5fd/mushroom.csv", "mushroom.csv")
    d = CSV.read("mushroom.csv", DataFrame)
    select!(d, [(atr => (x -> (;[Symbol(atr * "=" * string(value)) => d[!,atr] .== value for value in unique(d[!,atr])]...)) => AsTable)  for atr in names(d)])
end

function pretty(rules::Vector{Pair{Tuple{Vector{Symbol},Vector{Symbol}},Tuple{Float64,Float64}}})
    DataFrame([(antecedent=join(r[1][1],", "), consequent=join(r[1][2],", "), support=r[2][1], confidence=r[2][2])  for r in rules])
end

export apriori, dummy_dataset, dummy_dataset_biased, apriori_rule_gen, apriori_frequent_itemsets, merge_vectors, mushroom_dataset, pretty

end


text:
using TreatmentPanels
using Test
using DataFrames, Dates

@testset "Input tests" begin
    test_df = DataFrame(id = ["a", "a", "b", "b"], period = [1, 2, 1, 2], value = 1.0:4.0)

    # Column not present in data
    @test_throws ArgumentError BalancedPanel(test_df, "a" => 2; id_var = :wrong_var, t_var = :period, outcome_var = :value)
    @test_throws ArgumentError BalancedPanel(test_df, "a" => 2; id_var = :id, t_var = :wrong_var, outcome_var = :value)
    @test_throws "" BalancedPanel(test_df, "a" => 2; id_var = :id, t_var = :period, outcome_var = :wrong_var)

    # Treatment assignment has nonexistent units/time periods
    @test_throws "" BalancedPanel(test_df, "c" => 2; id_var = :id, t_var = :period, outcome_var = :value)
end

@testset "Single unit, continuous treatment" begin
    int_treatment = "Single Treated Unit" => 5
    date_treatment = "Single Treated Unit" => Date(2004)

    single_continuous_data = DataFrame(
        name = [fill("Single Treated Unit", 8); fill("Untreated Unit 1", 8); fill("Untreated Unit 2", 8)],
        period = repeat(1:8, 3),
        year_period = repeat(Date(2000):Year(1):Date(2007), 3),
        outcome = vec([parse(Int, "$(i)$(t)") for t ∈ 1:8, i ∈ 1:3]))

    # Treatment specified as a single pair
    sc_bp_int = BalancedPanel(single_continuous_data, int_treatment;
        id_var = :name, t_var = :period, outcome_var = :outcome)
    @test sc_bp_int isa BalancedPanel{SingleUnitTreatment{Continuous}}
    # Treatment specified as a length one vector of pairs

    sc_bp_int2 = BalancedPanel(single_continuous_data, [int_treatment];
        id_var = :name, t_var = :period, outcome_var = :outcome)
    @test sc_bp_int2 isa BalancedPanel{SingleUnitTreatment{Continuous}}

    # Year treatment
    sc_bp_year = BalancedPanel(single_continuous_data, date_treatment;
        id_var = :name, t_var = :year_period, outcome_var = :outcome)
    @test sc_bp_year isa BalancedPanel{SingleUnitTreatment{Continuous}}

    # Utility functions
    y₁₀, y₁₁, y₀₀, y₀₁ = decompose_y(sc_bp_year)
    @test y₁₀ == [11.0, 12.0, 13.0, 14.0]
    @test y₁₁ == [15.0, 16.0, 17.0, 18.0]
    @test y₀₀ == [21.0 22.0 23.0 24.0
                  31.0 32.0 33.0 34.0]
    @test y₀₁ == [25.0 26.0 27.0 28.0
                  35.0 36.0 37.0 38.0]

    @test length_T₀(sc_bp_year) == 4
    @test length_T₁(sc_bp_year) == 4
    @test treated_ids(sc_bp_year) == 1
    @test treated_labels(sc_bp_year) == "Single Treated Unit"
    @test first_treated_period_ids(sc_bp_year) == 5
    @test first_treated_period_ids(sc_bp_int) == 5
    @test first_treated_period_labels(sc_bp_year) == Date(2004)
end

@testset "Single unit, single time-limited treatment" begin
    treatment = "A" => Date(2000) => Date(2002)

    test_df = DataFrame(id = repeat(["A", "B", "C"], inner = 10),
        period = repeat(Date(2000):Year(1):Date(2009), 3),
        value = [collect(1:10); collect(2:11); collect(1.5:10.5)])

    bp = BalancedPanel(test_df, treatment;
        id_var = :id, t_var = :period, outcome_var = :value)

    @test bp isa BalancedPanel{SingleUnitTreatment{Discontinuous}}

    #!# TO DO - test utility functions for this case
end

@testset "Single unit, multiple time-limited treatments" begin
    #!# TO DO
    #treatment = "A" => [Date(2000) => Date(2002), Date(2005) => Date(2007)]

    #test_df = DataFrame(id = repeat(["A", "B", "C"], inner = 10),
    #    period = repeat(Date(2000):Year(1):Date(2009), 3),
    #    value = [collect(1:10); collect(2:11); collect(1.5:10.5)])

    #bp = BalancedPanel(test_df, treatment;
    #    id_var = :id, t_var = :period, outcome_var = :value)

    #@test bp isa BalancedPanel{SingleUnitTreatment, StartEndTreatment}
    @test true
end

@testset "Multiple units, continuous simultaneous treatment" begin
    treatment = ["A" => Date(2000), "B" => Date(2000)]

    test_df = DataFrame(id = repeat(["A", "B", "C"], inner = 10),
        period = repeat(Date(2000):Year(1):Date(2009), 3),
        value = [collect(1:10); collect(2:11); collect(1.5:10.5)])

    bp = BalancedPanel(test_df, treatment;
        id_var = :id, t_var = :period, outcome_var = :value)

    @test bp isa BalancedPanel{MultiUnitTreatment{Simultaneous{Continuous}}}

    #!# TO DO - test utility functions for this case
end

@testset "Multiple units, continuous staggered treatment" begin
    treatment = ["A" => Date(2000), "B" => Date(2001)]

    test_df = DataFrame(id = repeat(["A", "B", "C"], inner = 10),
        period = repeat(Date(2000):Year(1):Date(2009), 3),
        value = [collect(1:10); collect(2:11); collect(1.5:10.5)])

    bp = BalancedPanel(test_df, treatment;
        id_var = :id, t_var = :period, outcome_var = :value)

    @test bp isa BalancedPanel{MultiUnitTreatment{Staggered{Continuous}}}

    #!# TO DO - test utility functions for this case
end

#!# TO ADD - tests of plotting functionality