text:
mutable struct PSIContainer
    JuMPmodel::Union{Nothing, JuMP.AbstractModel}
    time_steps::UnitRange{Int}
    resolution::Dates.TimePeriod
    settings::PSISettings
    settings_copy::PSISettings
    variables::Dict{Symbol, AbstractArray}
    constraints::Dict{Symbol, AbstractArray}
    cost_function::JuMP.AbstractJuMPScalar
    expressions::Dict{Symbol, JuMP.Containers.DenseAxisArray}
    parameters::Union{Nothing, ParametersContainer}
    initial_conditions::InitialConditions
    pm::Union{Nothing, PM.AbstractPowerModel}

    function PSIContainer(
        sys::PSY.System,
        settings::PSISettings,
        jump_model::Union{Nothing, JuMP.AbstractModel},
    )
        PSY.check_forecast_consistency(sys)
        resolution = PSY.get_forecasts_resolution(sys)
        resolution = IS.time_period_conversion(resolution)
        new(
            jump_model,
            1:1,
            resolution,
            settings,
            copy_for_serialization(settings),
            Dict{Symbol, AbstractArray}(),
            Dict{Symbol, AbstractArray}(),
            zero(JuMP.GenericAffExpr{Float64, JuMP.VariableRef}),
            DenseAxisArrayContainer(),
            nothing,
            InitialConditions(use_parameters = get_use_parameters(settings)),
            nothing,
        )
    end
end

function PSIContainer(
    ::Type{T},
    sys::PSY.System,
    settings::PSISettings,
    jump_model::Union{Nothing, JuMP.AbstractModel},
) where {T <: PM.AbstractPowerModel}

    container = PSIContainer(sys, settings, jump_model)
    psi_container_init!(container, T, sys)
    return container
end

function _check_warm_start_support(JuMPmodel::JuMP.AbstractModel, warm_start_enabled::Bool)
    !warm_start_enabled && return warm_start_enabled
    solver_supports_warm_start =
        MOI.supports(JuMP.backend(JuMPmodel), MOI.VariablePrimalStart(), MOI.VariableIndex)
    if !solver_supports_warm_start
        solver_name = JuMP.solver_name(JuMPmodel)
        @warn("$(solver_name) does not support warm start")
    end
    return solver_supports_warm_start
end

function _make_jump_model!(psi_container::PSIContainer)
    settings = psi_container.settings
    parameters = get_use_parameters(settings)
    optimizer = get_optimizer(settings)
    if !isnothing(psi_container.JuMPmodel)
        if parameters
            if !haskey(psi_container.JuMPmodel.ext, :params)
                @info("Model doesn't have Parameters enabled. Parameters will be enabled")
                PJ.enable_parameters(psi_container.JuMPmodel)
                warm_start_enabled = get_warm_start(settings)
                solver_supports_warm_start =
                    _check_warm_start_support(psi_container.JuMPmodel, warm_start_enabled)
                set_warm_start!(settings, solver_supports_warm_start)
            end
        end
        return
    end
    @debug "Instantiating the JuMP model"
    if !isnothing(optimizer)
        JuMPmodel = JuMP.Model(optimizer)
        warm_start_enabled = get_warm_start(settings)
        solver_supports_warm_start =
            _check_warm_start_support(JuMPmodel, warm_start_enabled)
        set_warm_start!(settings, solver_supports_warm_start)
        parameters && PJ.enable_parameters(JuMPmodel)
        psi_container.JuMPmodel = JuMPmodel
    else
        @debug "The optimization model has no optimizer attached"
        JuMPmodel = JuMP.Model()
        parameters && PJ.enable_parameters(JuMPmodel)
        psi_container.JuMPmodel = JuMPmodel
    end
    return
end

function _make_container_array(parameters::Bool, ax...)
    if parameters
        return JuMP.Containers.DenseAxisArray{PGAE}(undef, ax...)
    else
        return JuMP.Containers.DenseAxisArray{GAE}(undef, ax...)
    end
    return
end

function _make_expressions_dict!(
    psi_container::PSIContainer,
    bus_numbers::Vector{Int},
    transmission::Type{S},
) where {S <: PM.AbstractPowerModel}
    settings = psi_container.settings
    parameters = get_use_parameters(settings)
    time_steps = 1:get_horizon(settings)
    psi_container.expressions = DenseAxisArrayContainer(
        :nodal_balance_active =>
            _make_container_array(parameters, bus_numbers, time_steps),
        :nodal_balance_reactive =>
            _make_container_array(parameters, bus_numbers, time_steps),
    )
    return
end

function _make_expressions_dict!(
    psi_container::PSIContainer,
    bus_numbers::Vector{Int},
    transmission::Type{S},
) where {S <: PM.AbstractActivePowerModel}
    settings = psi_container.settings
    parameters = get_use_parameters(settings)
    time_steps = 1:get_horizon(settings)
    psi_container.expressions = DenseAxisArrayContainer(
        :nodal_balance_active =>
            _make_container_array(parameters, bus_numbers, time_steps),
    )
    return
end

function psi_container_init!(
    psi_container::PSIContainer,
    ::Type{T},
    sys::PSY.System,
) where {T <: PM.AbstractPowerModel}
    PSY.set_units_base_system!(sys, "SYSTEM_BASE")
    # The order of operations matter
    settings = psi_container.settings
    _make_jump_model!(psi_container)
    @assert !isnothing(psi_container.JuMPmodel)
    make_parameters_container = get_use_parameters(settings)
    make_parameters_container && (psi_container.parameters = ParametersContainer())

    use_forecasts = get_use_forecast_data(settings)
    if make_parameters_container && !use_forecasts
        throw(IS.ConflictingInputsError("enabling parameters without forecasts is not supported"))
    end

    if get_initial_time(settings) == UNSET_INI_TIME
        set_initial_time!(settings, PSY.get_forecasts_initial_time(sys))
    end

    if get_horizon(settings) == UNSET_HORIZON
        set_horizon!(settings, PSY.get_forecasts_horizon(sys))
    end

    if use_forecasts
        total_number_of_devices = length(get_available_components(PSY.Device, sys))
        psi_container.time_steps = 1:get_horizon(settings)
        # The 10e6 limit is based on the sizes of the lp benchmark problems http://plato.asu.edu/ftp/lpcom.html
        # The maximum numbers of constraints and variables in the benchmark problems is 1,918,399 and 1,259,121,
        # respectively. See also https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2013/138847.pdf
        variable_count_estimate = length(psi_container.time_steps) * total_number_of_devices
        if variable_count_estimate > 10e6
            @warn("The estimated total number of variables that will be created in the model is $(variable_count_estimate). The total number of variables might be larger than 10e6 and could lead to large build or solve times.")
        end
    end

    bus_numbers = sort([PSY.get_number(b) for b in PSY.get_components(PSY.Bus, sys)])
    _make_expressions_dict!(psi_container, bus_numbers, T)
    return
end

function has_initial_conditions(psi_container::PSIContainer, key::ICKey)
    return has_initial_conditions(psi_container.initial_conditions, key)
end

function iterate_initial_conditions(psi_container::PSIContainer)
    return iterate_initial_conditions(psi_container.initial_conditions)
end

function get_initial_conditions(
    psi_container::PSIContainer,
    ::Type{T},
    ::Type{D},
) where {T <: InitialConditionType, D <: PSY.Device}
    return get_initial_conditions(psi_container, ICKey(T, D))
end

function get_initial_conditions(psi_container::PSIContainer, key::ICKey)
    return get_initial_conditions(psi_container.initial_conditions, key)
end

function set_initial_conditions!(psi_container::PSIContainer, key::ICKey, value)
    set_initial_conditions!(psi_container.initial_conditions, key, value)
end

_variable_type(cm::PSIContainer) = JuMP.variable_type(cm.JuMPmodel)
model_time_steps(psi_container::PSIContainer) = psi_container.time_steps
model_resolution(psi_container::PSIContainer) = psi_container.resolution
model_has_parameters(psi_container::PSIContainer) =
    get_use_parameters(psi_container.settings)
model_uses_forecasts(psi_container::PSIContainer) =
    get_use_forecast_data(psi_container.settings)
model_initial_time(psi_container::PSIContainer) = get_initial_time(psi_container.settings)
#Internal Variables, Constraints and Parameters accessors
get_variables(psi_container::PSIContainer) = psi_container.variables
get_constraints(psi_container::PSIContainer) = psi_container.constraints
get_parameters(psi_container::PSIContainer) = psi_container.parameters
get_expression(psi_container::PSIContainer, name::Symbol) = psi_container.expressions[name]
get_initial_conditions(psi_container::PSIContainer) = psi_container.initial_conditions
get_PTDF(psi_container::PSIContainer) = get_PTDF(psi_container.settings)
get_settings(psi_container::PSIContainer) = psi_container.settings
get_jump_model(psi_container::PSIContainer) = psi_container.JuMPmodel

function get_variable(
    psi_container::PSIContainer,
    var_type::AbstractString,
    ::Type{T},
) where {T <: PSY.Component}
    return get_variable(psi_container, make_variable_name(var_type, T))
end

function get_variable(
    psi_container::PSIContainer,
    ::Type{T},
    ::Type{U},
) where {T <: VariableType, U <: PSY.Component}
    return get_variable(psi_container, make_variable_name(T, U))
end

function get_variable(psi_container::PSIContainer, var_type::AbstractString)
    return get_variable(psi_container, make_variable_name(var_type))
end

function get_variable(psi_container::PSIContainer, update_ref::UpdateRef)
    return get_variable(psi_container, update_ref.access_ref)
end

function get_variable(psi_container::PSIContainer, name::Symbol)
    var = get(psi_container.variables, name, nothing)
    if isnothing(var)
        @error "$name is not stored" sort!(get_variable_names(psi_container))
        throw(IS.InvalidValue("variable $name is not stored"))
    end

    return var
end

function get_variable_names(psi_container::PSIContainer)
    return collect(keys(psi_container.variables))
end

function assign_variable!(
    psi_container::PSIContainer,
    variable_type::AbstractString,
    ::Type{T},
    value,
) where {T <: PSY.Component}
    assign_variable!(psi_container, make_variable_name(variable_type, T), value)
    return
end

function assign_variable!(psi_container::PSIContainer, variable_type::AbstractString, value)
    assign_variable!(psi_container, make_variable_name(variable_type), value)
    return
end

function assign_variable!(psi_container::PSIContainer, name::Symbol, value)
    @debug "assign_variable" name

    if haskey(psi_container.variables, name)
        @error "variable $name is already stored" sort!(get_variable_names(psi_container))
        throw(IS.InvalidValue("variable $name is already stored"))
    end

    psi_container.variables[name] = value
    return
end

function add_var_container!(
    psi_container::PSIContainer,
    var_name::Symbol,
    axs...;
    sparse = false,
)
    if sparse
        container = sparse_container_spec(psi_container.JuMPmodel, axs...)
    else
        container = container_spec(psi_container.JuMPmodel, axs...)
    end
    assign_variable!(psi_container, var_name, container)
    return container
end

function get_constraint(
    psi_container::PSIContainer,
    constraint_type::AbstractString,
    ::Type{T},
) where {T <: PSY.Component}
    return get_constraint(psi_container, make_constraint_name(constraint_type, T))
end

function get_constraint(psi_container::PSIContainer, constraint_type::AbstractString)
    return get_constraint(psi_container, make_constraint_name(constraint_type))
end

function get_constraint(psi_container::PSIContainer, name::Symbol)
    var = get(psi_container.constraints, name, nothing)
    if isnothing(var)
        @error "$name is not stored" sort!(get_constraint_names(psi_container))
        throw(IS.InvalidValue("constraint $name is not stored"))
    end

    return var
end

function get_constraint_names(psi_container::PSIContainer)
    return collect(keys(psi_container.constraints))
end

function assign_constraint!(
    psi_container::PSIContainer,
    constraint_type::AbstractString,
    ::Type{T},
    value,
) where {T <: PSY.Component}
    assign_constraint!(psi_container, make_constraint_name(constraint_type, T), value)
    return
end

function assign_constraint!(
    psi_container::PSIContainer,
    constraint_type::AbstractString,
    value,
)
    assign_constraint!(psi_container, make_constraint_name(constraint_type), value)
    return
end

function assign_constraint!(psi_container::PSIContainer, name::Symbol, value)
    @debug "set_constraint" name
    psi_container.constraints[name] = value
    return
end

function add_cons_container!(
    psi_container::PSIContainer,
    cons_name::Symbol,
    axs...;
    sparse = false,
)
    if sparse
        container = sparse_container_spec(psi_container.JuMPmodel, axs...)
    else
        container = JuMPConstraintArray(undef, axs...)
    end
    assign_constraint!(psi_container, cons_name, container)
    return container
end

function get_parameter_names(psi_container::PSIContainer)
    return collect(keys(psi_container.parameters))
end

function get_parameter_container(psi_container::PSIContainer, name::AbstractString)
    return get_parameter_container(psi_container, Symbol(name))
end

function get_parameter_container(psi_container::PSIContainer, name::Symbol)
    container = get(psi_container.parameters, name, nothing)
    if isnothing(container)
        @error "$name is not stored" sort!(get_parameter_names(psi_container))
        throw(IS.InvalidValue("parameter $name is not stored"))
    end
    return container
end

function get_parameter_container(
    psi_container::PSIContainer,
    name::Symbol,
    ::Type{T},
) where {T <: PSY.Component}
    return get_parameter_container(psi_container, encode_symbol(T, name))
end

function get_parameter_container(psi_container::PSIContainer, ref::UpdateRef)
    return get_parameter_container(psi_container, ref.access_ref)
end

function get_parameter_array(psi_container::PSIContainer, ref)
    return get_parameter_array(get_parameter_container(psi_container, ref))
end

function assign_parameter!(psi_container::PSIContainer, container::ParameterContainer)
    @debug "assign_parameter" container.update_ref
    name = container.update_ref.access_ref
    if name isa AbstractString
        name = Symbol(name)
    end

    if haskey(psi_container.parameters, name)
        @error "parameter $name is already stored" sort!(get_parameter_names(psi_container))
        throw(IS.InvalidValue("parameter $name is already stored"))
    end

    psi_container.parameters[name] = container
    return
end

function add_param_container!(
    psi_container::PSIContainer,
    param_reference::UpdateRef,
    axs...,
)
    container = ParameterContainer(
        param_reference,
        JuMP.Containers.DenseAxisArray{PJ.ParameterRef}(undef, axs...),
        JuMP.Containers.DenseAxisArray{Float64}(undef, axs...),
    )
    assign_parameter!(psi_container, container)
    return container
end

function iterate_parameter_containers(psi_container::PSIContainer)
    Channel() do channel
        for container in values(psi_container.parameters)
            put!(channel, container)
        end
    end
end

function get_parameters_value(psi_container::PSIContainer)
    # TODO: Still not obvious implementation since it needs to get the multipliers from
    # the system
    params_dict = Dict{Symbol, DataFrames.DataFrame}()
    parameters = get_parameters(psi_container)
    (isnothing(parameters) || isempty(parameters)) && return params_dict
    for (k, v) in parameters
        !isa(v.update_ref, UpdateRef{<:PSY.Component}) && continue
        params_key_tuple = decode_symbol(k)
        params_dict_key = Symbol(params_key_tuple[1], "_", params_key_tuple[3])
        param_array = axis_array_to_dataframe(get_parameter_array(v))
        multiplier_array = axis_array_to_dataframe(get_multiplier_array(v))
        params_dict[params_dict_key] = param_array .* multiplier_array
    end
    return params_dict
end

function assign_expression!(psi_container::PSIContainer, name::Symbol, value)
    @debug "set_expression" name
    psi_container.expressions[name] = value
    return
end

function add_expression_container!(psi_container::PSIContainer, exp_name::Symbol, axs...)
    container = JuMP.Containers.DenseAxisArray{JuMP.GenericAffExpr}(undef, axs...)
    assign_expression!(psi_container, exp_name, container)
    return container
end

function is_milp(container::PSIContainer)
    type_of_optimizer = typeof(container.JuMPmodel.moi_backend.optimizer.model)
    supports_milp = hasfield(type_of_optimizer, :last_solved_by_mip)
    !supports_milp && return false
    return container.JuMPmodel.moi_backend.optimizer.model.last_solved_by_mip
end

function _export_optimizer_log(
    optimizer_log::Dict{Symbol, Any},
    psi_container::PSIContainer,
    path::String,
)
    optimizer_log[:obj_value] = JuMP.objective_value(psi_container.JuMPmodel)
    optimizer_log[:termination_status] =
        Int(JuMP.termination_status(psi_container.JuMPmodel))
    optimizer_log[:primal_status] = Int(JuMP.primal_status(psi_container.JuMPmodel))
    optimizer_log[:dual_status] = Int(JuMP.dual_status(psi_container.JuMPmodel))
    try
        optimizer_log[:solve_time] = MOI.get(psi_container.JuMPmodel, MOI.SolveTime())
    catch
        @warn("SolveTime() property not supported by the Solver")
        optimizer_log[:solve_time] = NaN # "Not Supported by solver"
    end
    write_optimizer_log(optimizer_log, path)
    return
end

""" Exports the OpModel JuMP object in MathOptFormat"""
function _write_psi_container(psi_container::PSIContainer, save_path::String)
    MOF_model = MOPFM(format = MOI.FileFormats.FORMAT_MOF)
    MOI.copy_to(MOF_model, JuMP.backend(psi_container.JuMPmodel))
    MOI.write_to_file(MOF_model, save_path)
    return
end

function get_dual_values(psi_container::PSIContainer)
    cons = get_constraint_duals(psi_container.settings)
    return get_dual_values(psi_container, cons)
end

function get_dual_values(op::PSIContainer, cons::Vector{Symbol})
    results_dict = Dict{Symbol, DataFrames.DataFrame}()
    isempty(cons) && return results_dict
    for c in cons
        v = get_constraint(op, c)
        results_dict[c] = axis_array_to_dataframe(v)
    end
    return results_dict
end

function add_to_setting_ext!(psi_container::PSIContainer, key::String, value)
    settings = get_settings(psi_container)
    push!(get_ext(settings), key => value)
    @debug "Add to settings ext" key value
    return
end


text:
include("overapprox_1d.jl")
include("autoline.jl")
include("overt_utils.jl")
include("definitions.jl")
using SymEngine
using Plots
plotly()
using PGFPlots

# TODO: put params into a mutable struct
plotflag = false
plottype = "pgf"
DELETE_MUL_BY_ZERO = true
DELETE_DEAD_RELUS = true

"""
    overapprox(expr,
                  range_dict::Dict{Symbol, Array{T, 1}} where {T <: Real};
                  N::Integer=2,
                  ϵ::Real=1e-2)
    -> OverApproximation
Overapproximate an n-dimensional function using a relational abstraction.

N    is the number of linear segments used per region of constant convexity. To allow OVERT to automatically choose this number, set N to -1. It will add segments until the maximum deviation between the function and the bound is within 2%.

ϵ    is the RELATIVE "air gap" added to the bounds so that the bounds don't touch the function at the point of closest approach. If the function bound lies in the range [a,b] a gap of size ϵ*(b-a) will be added to the y values of all points of the upper bound and subtracted from all points of the lower bound. The default value is 1% or ϵ = 0.01
"""
function overapprox(expr,
                       range_dict::Dict{Symbol, Array{T, 1}} where {T <: Real};
                       N::Integer=2,
                       ϵ::Real=1e-2)
    println("Using N=$N, ϵ=$ϵ")
    bound = OverApproximation()
    range_dict = floatize(range_dict)
    bound.ranges = range_dict
    bound.N = N
    bound.ϵ = ϵ
    return overapprox(expr, bound)
end

"""
    overapprox(expr, bound::OverApproximation) -> bound::OverApproximation

Helper function. Bounds the expression _expr_ in a recursive manner. bound is used as a container
to collect the relations describing the bound along the way.
"""
function overapprox(expr,
                       bound::OverApproximation)
    # all operations should have at most two arguments.
    expr = reduce_args_to_2(expr)
    @debug(expr)

    # base cases
    if expr isa Symbol
        @debug("symbol base case")
        bound.output = expr
        bound.output_range = bound.ranges[expr]
        return bound
    elseif is_number(expr)
        @debug("is number base case")
        bound.output = eval(expr)
        return bound
        ### This code below preserves an exact numeric expression but adds variables
        # bound.output = add_var(bound)
        # c  = eval(expr)
        # bound.output_range = [c,c]
        # # use :($newvar = $expr) # for use with infinite precision solvers
        # push!(bound.approx_eq, :($(bound.output) == $c))
        # bound.fun_eq[bound.output] = c
        # bound.ranges[bound.output] = bound.output_range
        # push!(bound.consts, bound.output)
        # return bound
    elseif is_affine(expr)
        @debug("$expr is affine base case")
        bound.output = add_var(bound)
        bound.output_range = [find_affine_range(expr, bound.ranges)...]
        push!(bound.approx_eq, :($(bound.output) == $expr))
        bound.fun_eq[bound.output] = expr
        @debug "turning $expr into $(bound.output) = $expr"
        bound.ranges[bound.output] = bound.output_range
        return bound
    # elseif is_1d(expr)
    #     """
    #     Currently, this case would resort to relying on fzero to find roots of d2f which
    #     is not sound as this is a numerical procedure. Thus it has been removed.
    #     TODO in the future is to implement a sound version of this for arbitrary
    #     1D functions using symbolic differentiation from Calculus.jl:differentiate to find the
    #     d2f function and then use the package IntervalRootFinding to find all roots in the
    #     interval in question in a GUARANTEED way!
    #     Allowing the usage of this case (is_1d) can speed up verification with OVERT for some functions by
    #     reducing the number of approximations needed and thus reducing the number of relus
    #     used. E.g. sin(x) + log(x^2) is a 1D function but without using the is_1d case, this
    #     would incur 3 separate approximations: x^2, log(v_3) and sin(x) . If this case can be used,
    #     it will recognize that sin(x) + log(x^2) is a 1D function and only construct a single
    #     OVERT approximation.
    #     """
    #     @debug("is 1d base case")
    #     sym = find_variables(expr)
    #     arg = sym[1]
    #     bound = bound_1arg_function(expr, arg, bound)
    #     return bound
    # recursive cases
    else # is a non-affine function
        if is_unary(expr)
            @debug("is unary")
            # unary expression with two args, like sin(x), cos(x), ...
            f = expr.args[1]
            xexpr = expr.args[2]
            # recurse on argument
            bound = overapprox(xexpr, bound)
            if f == :- # handle unary minus
                @debug "let affine handle unary minus"
                new_expr = :(-$(bound.output))
                return overapprox(new_expr, bound)
            elseif f == :relu
                @debug "PWL case"
                bound = handle_PWL(f, bound.output, bound)
                return bound
            end
            # handle outer function f
            bound = bound_unary_function(f, bound)
            return bound
        else # is binary
            # pull out the big guns
            @assert is_binary(expr)
            @debug("is binary")
            f = expr.args[1]
            xexpr = expr.args[2]
            yexpr = expr.args[3]
            # recurse on aguments
            bound = overapprox(xexpr, bound)
            x = bound.output
            #
            bound = overapprox(yexpr, bound)
            y = bound.output
            # handle outer function f
            bound = bound_binary_functions(f, x, y, bound)
            return bound
        end
    end
end

function bound_binary_functions(f, x, y, bound) # TODO: should only be for when both args are vars or 6/x
    @debug "bound binary functions: f x y:" f x y
    mul_two_vars = !is_number(x) && !(is_number(y)) && f  == :*
    divide_by_var = !is_number(y) && f == :/
    divide_by_const = is_number(y) && f ==:/
    if mul_two_vars
        # both args contain variables
        @debug("bounding * btw two variables")
        # first expand multiplication expr
        expanded_mul_expr, bound = expand_multiplication(x, y, bound)
        # and then bound that expression wrt x2, y2
        bound = overapprox(expanded_mul_expr, bound)
        return bound
    elseif divide_by_var
        # handles x/y and c/y
        # x/c should be transformed into (1./c)*x by simplify
        @debug("division x/y or c/y")
        # range of y should not include 0
        if bound.ranges[y][1]*bound.ranges[y][2] <= 0
            error("range of y should not include 0")
        end
        if is_number(x) # c/y
            @debug("c/y")
            f_new = :($x/$y)
            bound = bound_1arg_function(f_new, y, bound)
            return bound
        elseif x isa Symbol # x/y -> x*(1/y)
            @debug "converting: $x / $y"
            inv_denom = :(1. / $y)
            f_new = :($x * $inv_denom)
            @debug "into: $f_new"
            return overapprox(f_new, bound)
        end
    elseif divide_by_const
        new_expr = rewrite_division_by_const(:($f($x,$y)))
        return overapprox(new_expr, bound)
    # The following operation is not yet fully tested.
    elseif f == :^ # f(x,y)   x^2
        if is_number(y) # this is f(x)^a, where a is a constant.
            @debug "handling f(x)^a"
            f_new = :($x^$y)
            bound = bound_1arg_function(f_new, x, bound)
            return bound
        elseif is_number(x)  # this is a^f(x), where a is a constant.
            @debug "handling  a^f(x)"
            f_new = :($x^$y)
            bound = bound_1arg_function(f_new, y, bound)
            return bound
        else  # this is f(x)^g(y)
            error("f(x)^g(y) is not implemented yet.")
            return OverApproximation() # type stability?
        end
    elseif is_affine(:($f($x,$y)))
        new_expr = :($f($x,$y))
        @debug "Recursing to let affine handle: " new_expr
        return overapprox(new_expr, bound)
    elseif f ∈ [:min, :max, :relu]
        @debug "PWL case, f=" f
        bound = handle_PWL(f, x, y, bound)
        return bound
    else
        error(" operation $f with operands $x and $y is not implemented")
        return OverApproximation() # for type stability, maybe?
    end
end

function handle_PWL(f, x, y, bound::OverApproximation)
    # two arg case
    f_range = find_range(:($f($x, $y)), bound.ranges)
    bound.output_range = f_range
    if f_range == [0.0, 0.0] && DELETE_DEAD_RELUS # dead relu
        println("Deleting dead relu. f_range is $(f_range). ")
        bound.output = 0.0
    else # live relu or max or something
        bound.output = add_var(bound)
        expr = :($(bound.output) == $f($x, $y))
        push!(bound.approx_eq, expr)
        bound.fun_eq[bound.output] = :($f($x, $y))
        @debug "turning $f($x,$y) into $(bound.output) = $(bound.fun_eq[bound.output])"
        bound.ranges[bound.output] = bound.output_range
    end
    return bound
end

function handle_PWL(f, x, bound::OverApproximation)
    # one arg case
    f_range = find_range(:($f($x)), bound.ranges)
    bound.output_range = f_range
    if f_range == [0.0, 0.0] && DELETE_DEAD_RELUS # dead relu
        println("Deleting dead relu. f_range is $(f_range). ")
        bound.output = 0.0
    else # live relu or something
        bound.output = add_var(bound)
        expr = :($(bound.output) == $f($x))
        push!(bound.approx_eq, expr)
        bound.fun_eq[bound.output] = :($f($x))
        @debug "turning $f($x) into $(bound.output) = $(bound.fun_eq[bound.output])"
        bound.ranges[bound.output] = bound.output_range
    end
    return bound
end

function expand_multiplication(x, y, bound; ξ=0.1)
    return expand_multiplication_with_scaling(x, y, bound; ξ=ξ)
end

function expand_multiplication_with_scaling(x, y, bound; ξ=0.1)
    """
        expand_multiplication(x, y, bound; ξ=1.0)
    Re write multiplication e.g. x*y using exp(log()) and affine expressions
    e.g. x*y, x ∈ [a,b] ∧ y ∈ [c,d], ξ>0
         x2 = (x - a)/(b - a) + ξ  , x2 ∈ [ξ, 1 + ξ] aka x2 > 0   (recall, b > a)
            x = (b - a)(x2 - ξ) + a
         y2 = (y - c)/(d - c) + ξ , y2 ∈ [ξ, 1 + ξ] aka y2 > 0   (recall, d > c)
            y = (d - c)(y2 - ξ) + c
        x*y = ((b - a)(x2 - ξ) + a )*((d - c)(y2 - ξ) + c)
            = (b - a)*(d - c)*x2*y2  + (b - a)(c - ξ(d - c))*x2 + (d - c)(a - ξ(b - a))*y2 + (c - ξ(d - c))(a - ξ(b - a))
            = (b - a)*(d - c)*exp(log(x2*y2)) + (b - a)(c - ξ(d - c))*x2 + (d - c)(a - ξ(b - a))*y2 + (c - ξ(d - c))(a - ξ(b - a))
            = (b - a)*(d - c)*exp(log(x2) + log(y2)) + (b - a)(c - ξ(d - c))*x2 + (d - c)(a - ξ(b - a))*y2 + (c - ξ(d - c))(a - ξ(b - a))
        In this final form, everything is decomposed into unary functions, +/-, and affine functions!
    """

    a,b = bound.ranges[x]
    c,d = bound.ranges[y]
    @debug "mult ranges: [$a, $b], [$c, $d]"

    if ((a==0 && b==0) || (c==0 && d==0)) && DELETE_MUL_BY_ZERO
        println("Multiplication by zero.")
        # NOTE: in the future... might want to perserve e.g. v_3 = 0*v_2 and set e.g. v_3 as the output with range 0.0 rather than setting the output to 0.0 directly ¯\_(ツ)_/¯ not sure
        bound.output = 0.0
        return 0.0, bound
    else
        @assert(b >= a)
        @assert(d >= c)
        b_minus_a = b - a
        d_minus_c = d - c
        x2 = add_var(bound)
        y2 = add_var(bound)
        push!(bound.approx_eq, :($x2 == ($x - $a)/$b_minus_a + $ξ))
        push!(bound.approx_eq, :($y2 == ($y - $c)/$d_minus_c + $ξ))
        @debug("Expanding multiplication")
        bound.fun_eq[x2] = :(($x - $a)/$b_minus_a + $ξ)
        bound.fun_eq[y2] = :(($y - $c)/$d_minus_c + $ξ)
        @debug "$x2 = $(bound.fun_eq[x2])"
        @debug "$y2 = $(bound.fun_eq[y2])"
        bound.ranges[x2] = [ξ, 1. + ξ]
        bound.ranges[y2] = [ξ, 1. + ξ]

        b_minus_a_times_d_minus_c = (b - a)*(d - c)
        mult_expr = :($b_minus_a_times_d_minus_c*exp(log($x2) + log($y2)))
        x_coeff = (b - a)*(c - ξ*(d - c))
        xterm = :($x_coeff*$x2)
        y_coeff = (d - c)*(a - ξ*(b - a))
        yterm = :($y_coeff*$y2)
        constant_term = (c - ξ*(d - c))*(a - ξ*(b - a))

        expr = :($mult_expr + $xterm + $yterm + $constant_term)
        @debug "expanded multiplication expr is: $expr"
        return expr, bound
    end
end

function apply_fx(f, a)
    substitute!(f, :x, a)
end

"""
    bound_1arg_function(e::Expr, x::Symbol, bound::OverApproximation; plotflag=true)
Bound one argument functions like x^2.
Create upper and lower bounds of function f(x)
"""
function bound_1arg_function(e::Expr, arg::Symbol, bound::OverApproximation; plotflag=plotflag)
    @debug "bound effectively unary" e arg bound
    fun = SymEngine.lambdify(e, [arg])
    lb, ub, npoint = bound.ranges[arg][1], bound.ranges[arg][2], bound.N
    d2f_zeros, convex = get_regions_1arg(e, arg, lb, ub)
    return bound_unary_function(fun, e, arg, lb, ub, npoint, bound, plotflag=plotflag, d2f_zeros=d2f_zeros, convex=convex)
end

"""
    bound_unary_function(f::Symbol, x_bound::OverflowError; plotflag=true)
Bound one argument unary functions like sin(x). Create upper and lower bounds of function f(x)
"""
function bound_unary_function(f::Symbol, x_bound::OverApproximation; plotflag=plotflag)
    @debug "bound true unary" f x_bound.output
    fun = eval(:($f))
    lb, ub, npoint = x_bound.output_range[1], x_bound.output_range[2], x_bound.N
    f_x_expr = :($(f)($(x_bound.output)))
    d2f_zeros, convex = get_regions_unary(f, lb, ub)
    return bound_unary_function(fun, f_x_expr, x_bound.output, lb, ub, npoint, x_bound, plotflag=plotflag, d2f_zeros=d2f_zeros, convex=convex)
end

"""
    bound_unary_function(f::Function, lb, ub, npoint, bound)
Bound one argument functions like sin(x) or x -> x^2 or x -> 1/x. Create upper and lower bounds of function f(x)
"""
function bound_unary_function(fun::Function, f_x_expr, x, lb, ub, npoint, bound; plotflag=plotflag, d2f_zeros=nothing, convex=nothing)
    UBpoints, UBfunc_sym = find_UB(fun, lb, ub, npoint; lb=false, plot=plotflag, ϵ= bound.ϵ, d2f_zeros=d2f_zeros, convex=convex)
    fUBrange = [find_1d_range(UBpoints)...]
    LBpoints, LBfunc_sym = find_UB(fun, lb, ub, npoint; lb=true, plot=plotflag, ϵ= -bound.ϵ, d2f_zeros=d2f_zeros, convex=convex)
    fLBrange = [find_1d_range(LBpoints)...]

    # plot after adding air gap and after turning into closed form expression
    @debug "bounding true unary..."
    if plotflag
        p = Plots.plot(0,0)
        global NPLOTS
        NPLOTS += 1
        if plottype != "pgf"
            p = plot(range(lb, ub, length=100), fun, label="function", color="red")
            plot!(p, [p[1] for p in LBpoints], [p[2] for p in LBpoints],  color="purple", marker=:o, markersize=1, label="lower bound")
            plot!(p, [p[1] for p in UBpoints], [p[2] for p in UBpoints], color="blue", marker=:diamond, markersize=1,  label="upper bound", legend=:right, title="Function and bounds")
            # display(p)
            savefig(p, "plots/bound_"*string(NPLOTS)*".html")
        else # plottype == pgf
            println("Saving PGF plot")
            if f_x_expr.args[1] ∈ [:/, :^]
                # use whole expression in title
                fun_string = "\$"*string(f_x_expr)*"\$"
            else
                fun_string = "\$"*string(fun)*"\$"
            end
            println("funstring = $(fun_string)")
            x_plot_points = range(lb, ub, length=100)
            f_x_plot_points = fun.(x_plot_points)
            fig = PGFPlots.Axis(style="width=10cm, height=10cm", ylabel="\$f(x)\$", xlabel="x", title=fun_string)
            push!(fig, PGFPlots.Plots.Linear(x_plot_points, f_x_plot_points, legendentry=fun_string, style="solid, black, mark=none"))
            push!(fig, PGFPlots.Plots.Linear([p[1] for p in LBpoints], [p[2] for p in LBpoints], legendentry="lower bound", style="solid, purple, mark=none"))
            push!(fig, PGFPlots.Plots.Linear([p[1] for p in UBpoints], [p[2] for p in UBpoints], legendentry="upper bound", style="solid, blue, mark=none"))
            fig.legendStyle =  "at={(1.05,1.0)}, anchor=north west"
            PGFPlots.save("plots/bound_"*string(NPLOTS)*".tex", fig)
            #PGFPlots.save("plots/bound_"*string(NPLOTS)*".pdf", fig)
        end

    end
    @debug "LB function is: $LBfunc_sym"
    @debug "UB function is: $UBfunc_sym"

    ## create new vars for these expr, equate to exprs, and add them to equality list
    # e.g. y = fUB(x), z = fLB(x)
    UBvar = add_var(bound)
    push!(bound.approx_eq, :($UBvar == $(apply_fx(UBfunc_sym, x))))
    bound.ranges[UBvar] = fUBrange
    LBvar = add_var(bound)
    push!(bound.approx_eq, :($LBvar == $(apply_fx(LBfunc_sym, x))))
    bound.ranges[LBvar] = fLBrange
    OArange = [min(fLBrange...,fUBrange...), max(fLBrange..., fUBrange...)]
    ## create variable representing overapprox and add inequality to ineq list
    OAvar = add_var(bound)
    push!(bound.approx_ineq, (:($LBvar ≦ $OAvar)) )
    push!(bound.approx_ineq, (:($OAvar ≦ $UBvar)) )
    bound.fun_eq[OAvar] = f_x_expr
    bound.ranges[OAvar] = OArange
    bound.output = OAvar
    bound.output_range = OArange
    return bound
end

function floatize(dict::Dict{Symbol, Array{T, 1}} where {T <: Real})
    # values of dict must be convertible to float
    newdict = Dict{Symbol,Array{Float64,1}}()
    for k=keys(dict)
        newdict[k] = float(dict[k])
    end
    return newdict::Dict{Symbol, Array{T, 1}} where {T <: Real}
end

ext: jl


text:
@enum Polarity::Bool DATA ANTIDATA

@enum WaiterState begin
    WAITER_INIT
    WAITER_WAITING
    WAITER_NOTIFYING
    WAITER_SATISFIED
    WAITER_FINISHED
    WAITER_TRASHED
    WAITER_INTERRUPTED
end

mutable struct Waiter{T}
    # _tag_pad::PadAfter32
    task::Union{Task,Nothing}
    value::Union{Nothing,T}
    @atomic state::WaiterState
    # _state_pad::PadAfter64  # TODO: smaller pad
end

Waiter{T}(task::Task = current_task()) where {T} = Waiter{T}(
    # _PadAfter32(),
    task,
    nothing,
    WAITER_INIT,
    # _PadAfter64(),
)

function reinit!(w::Waiter)
    if assertion_enabled()
        state = @atomic w.state
        @assert state in (WAITER_TRASHED, WAITER_FINISHED)
        @assert w.task === nothing
        @assert w.value === nothing
    end
    w.task = current_task()
    w.value === nothing
    @atomic :monotonic w.state = WAITER_INIT
    return w
end

function set_interrupted!(w::Waiter, nspins::Integer = 100)
    old, ok = @atomicreplace w.state WAITER_WAITING => WAITER_INTERRUPTED
    for _ in 1:nspins
        ok && break
        old, ok = @atomicreplace w.state old => WAITER_INTERRUPTED
        GC.safepoint()
    end
    if !ok
        @atomic w.state = WAITER_INTERRUPTED
    end
    if !ok || old != WAITER_WAITING
        # TODO: is `@warn` better?
        @debug "an item lost" ok old
    end
    return old, ok
end

function Base.fetch(w::Waiter{T}; nspins::Integer = 0) where {T}
    for _ in 1:nspins
        if (@atomic w.state) == WAITER_SATISFIED
            @goto satisifed
        end
        GC.safepoint()
    end

    _, success = @atomicreplace w.state WAITER_INIT => WAITER_WAITING
    if success
        try
            wait()
        catch
            set_interrupted!(w)
            rethrow()
        end
    end

    @label satisifed
    let old = @atomic w.state
        @assert old in (WAITER_SATISFIED, WAITER_NOTIFYING)
        old, success = @atomicreplace w.state old => WAITER_FINISHED
        if !success
            error("more than two duals detected: old = $old")
        end
    end

    x = w.value
    # Help GC:
    w.task = nothing
    w.value = nothing
    return x::T
end

function tryput!(w::Waiter{T}, x::T) where {T}
    # Before the CAS, `put!` owns the field `.value`:
    w.value = x

    old = WAITER_WAITING  # initial guess
    while true
        if old == WAITER_WAITING
            new = WAITER_NOTIFYING
        elseif old == WAITER_INTERRUPTED
            return false
        else
            @assert old == WAITER_INIT "required old == WAITER_INIT but old = $old"
            new = WAITER_SATISFIED
        end
        old, success = @atomicreplace w.state old => new
        if success
            if old == WAITER_WAITING
                schedule(w.task::Task)
            end
            return true
        end
    end
end

function trash!(w::Waiter)
    if assertion_enabled()
        old, success = @atomicreplace w.state WAITER_INIT => WAITER_TRASHED
        if !success
            error("unshared waiter mutaed: state: $old")
        end
    else
        @atomic :monotonic w.state = WAITER_TRASHED
    end
    w.task = nothing
end

struct MPCRQSlot{S}
    index_safe_polarity::UInt32
    storage::S  # sizeof(S) ≤ 4
end

const MPCRQSLOT_MAX_INDEX = Int32(typemax(UInt32) >> 2)

@inline function MPCRQSlot(; index::Integer, safe::Bool, polarity::Polarity, storage)
    # index ≤ typemax(UInt32) >> 2 || error("index $index too large")
    index_safe_polarity = (index % UInt32) << 2
    index_safe_polarity |= (safe % UInt32) << 1
    index_safe_polarity |= UInt32(polarity)
    return MPCRQSlot(index_safe_polarity, storage)
end

@inline MPCRQSlot{T}(;
    index::Integer,
    safe::Bool,
    polarity::Polarity,
    storage::T = zero(T),
) where {T<:Number} = MPCRQSlot(; index, safe, polarity, storage)

@inline function Base.getproperty(slot::MPCRQSlot, name::Symbol)
    index_safe_polarity = getfield(slot, :index_safe_polarity)
    if name === :index
        return (index_safe_polarity >> 2) % Int32
    elseif name === :safe
        return ((index_safe_polarity >> 1) & 0x01) % Bool
    elseif name === :polarity
        return Polarity((index_safe_polarity & 0x01) % Bool)
    end
    return getfield(slot, name)
end

@enum MPCRQResult MPCRQ_CLOSED MPCRQ_ENQUEUED

# See: IndirectConcurrentRingQueueNode
mutable struct IndirectMultiPolarityConcurrentRingQueueNode{
    T,
    Data<:Union{Nothing,Vector{Vector{T}}},
}
    _tag_pad::PadAfter32
    @atomic data_idx::Int32
    _data_idx_pad::PadAfter32
    @atomic antidata_idx::Int32
    _antidata_idx_pad::PadAfter32
    @atomic closed_info::Int32
    _closed_info_pad::PadAfter32
    @atomic next::Union{Nothing,IndirectMultiPolarityConcurrentRingQueueNode{T}}
    _next_pad::PadAfter64
    ring::typeof(cacheline_padded_vector(MPCRQSlot{UInt32}, 1))
    log2len::Int
    # buffers::Vector{Vector{Some{T},Waiter{T}}}
    data::Data
    waiters::Vector{Vector{Waiter{T}}}
end

inlinedata(::IndirectMultiPolarityConcurrentRingQueueNode{<:Any,Data}) where {Data} =
    Data === Nothing

# TODO: more types
const SmallSigned = Union{Int8,Int16}
const SmallUnsigned = Union{Bool,UInt8,UInt16}
const NullableInline32 = Union{SmallSigned,SmallUnsigned}

function IndirectMultiPolarityConcurrentRingQueueNode{T}(log2len::Int) where {T}
    len = 2^log2len
    if T <: NullableInline32 && isconcretetype(T)
        data = nothing
    else
        data = [Vector{T}(undef, len) for _ in 1:Threads.nthreads()]
    end
    waiters = [Vector{Waiter{T}}(undef, len) for _ in 1:Threads.nthreads()]
    ring = init_ring!(cacheline_padded_vector(MPCRQSlot{UInt32}, len))
    data_idx = Int32(1)
    antidata_idx = Int32(1)
    closed_info = Int32(1)
    return IndirectMultiPolarityConcurrentRingQueueNode(
        _PadAfter32(),
        data_idx,
        _PadAfter32(),
        antidata_idx,
        _PadAfter32(),
        closed_info,
        _PadAfter32(),
        nothing,
        _PadAfter64(),
        ring,
        log2len,
        data,
        waiters,
    )
end
# It looks like some non-neglegible amount of time is spent in
# `Vector{Waiter{T}}(undef, _)`, presumably zeroing out the elements.  A tricky
# part of this is detecting that the CRQ is no longer used. Since the CRQ may be
# closed just after the FAI and before the CAS on the slot is complete, (naively
# thinking) linear scan the `.ring` array is required for checking that the CRQ
# is not longer used (without adding other mechanisms like ref conting). So,
# ATM, only the newcrq that lost CAS is cached (which still helps the
# performance).

function init_ring!(ring)
    for index in eachindex(ring)
        ring[index] = MPCRQSlot{UInt32}(; index, safe = true, polarity = DATA)
    end
    return ring
end

function Base.empty!(crq::IndirectMultiPolarityConcurrentRingQueueNode)
    @atomic crq.data_idx = Int32(1)
    @atomic crq.antidata_idx = Int32(1)
    @atomic crq.closed_info = Int32(1)
    @atomic crq.next = nothing
    init_ring!(crq.ring)
    # for waiters in crq.waiters
    #     n = length(waiters)
    #     resize!(empty!(waiters), n)
    # end
    return crq
end

Base.similar(crq::IndirectMultiPolarityConcurrentRingQueueNode{T}) where {T} =
    IndirectMultiPolarityConcurrentRingQueueNode{T}(crq.log2len)

Base.eltype(::Type{<:IndirectMultiPolarityConcurrentRingQueueNode{T}}) where {T} = T

#=
function Base.isempty(crq::IndirectMultiPolarityConcurrentRingQueueNode)
    closed_info = @atomic crq.closed_info
    closed_info_idx = extract_index(closed_info)
    isclosed = closed_info < 0
    isclosed || return false
    data_idx = extract_index(@atomic crq.data_idx)
    data_idx ≥ closed_info_idx || return false
    antidata_idx = extract_index(@atomic crq.antidata_idx)
    antidata_idx ≥ closed_info_idx || return false
    # TODO: relaxed loads plus fence
    max_idx = max(data_idx, antidata_idx)
    for slot in crq.ring
        (; index, safe) = slot
        safe && index ≤ max_idx && return false
    end
    return true
end
=#

@inline function nonnegative(x)
    if assertion_enabled()
        @assert x ≥ 0
    end
    return x
end

@inline function negative(x)
    if assertion_enabled()
        @assert x < 0
    end
    return x
end

@inline set_closing(i) = negative(ifelse(i < zero(i), i, typemin(i) + i))
@inline extract_index(i) = nonnegative(ifelse(i ≥ zero(i), i, i - typemin(i)))

@inline _mod1(i::Int32, crq) = _mod1log2(i, crq.log2len)
@inline _mod1log2(i::Int32, log2x) =
    (UInt32(i - 1) & (typemax(UInt32) >> (32 - log2x))) % Int32 + Int32(1)

_embed(::Type{UInt32}, x::SmallSigned) = _embed(UInt32, unsigned(x))
_embed(::Type{UInt32}, x::SmallUnsigned) = (x % UInt32) + UInt32(1)
_extract(::Type{T}, x::UInt32) where {T<:SmallSigned} = signed(_extract(unsigned(T), x))::T
_extract(::Type{T}, x::UInt32) where {T<:SmallUnsigned} = (x - UInt32(1)) % T

"""
    denqueue!(crq::IMPCRQ{T}, x::T) -> MPCRQ_CLOSED or MPCRQ_ENQUEUED or y::Waiter
    denqueue!(crq::IMPCRQ{T}, x::Waiter{T}) -> MPCRQ_CLOSED or MPCRQ_ENQUEUED or y::Some{T}
    where MPCRQ = IndirectMultiPolarityConcurrentRingQueueNode
"""
function denqueue!(
    crq::IndirectMultiPolarityConcurrentRingQueueNode{T},
    x::Union{T,Waiter{T}},
) where {T}
    nthreads = Int32(Threads.nthreads())

    local polarity::Polarity = x isa Waiter{T} ? ANTIDATA : DATA
    R = 2^crq.log2len
    starvation_counter = 0
    while true
        if x isa Waiter{T}
            # p, _ = modifyfield!(crq, :antidata_idx, +, Int32(1), :sequentially_consistent)
            p = atomic_modifyfield!(crq, Val(:antidata_idx), +, true)::Int32
        else
            # p, _ = modifyfield!(crq, :data_idx, +, Int32(1), :sequentially_consistent)
            p = atomic_modifyfield!(crq, Val(:data_idx), +, true)::Int32
        end
        p_closing = p < 0
        p_idx = extract_index(p)
        if p_closing
            if discovered_closing(crq) ≤ p_idx
                return MPCRQ_CLOSED
            end
        elseif p_idx ≥ MPCRQSLOT_MAX_INDEX - nthreads
            discovered_closing(crq)
            error("NOT WELL TESTED")
            return MPCRQ_CLOSED
        end
        itemindex = _mod1(p_idx, crq)
        slotptr = pointer(crq.ring, itemindex)
        # ncas = 0
        while true
            # ncas += 1
            # if ncas > 10000
            #     let msg = "too many retires tid=$(Threads.threadid())"
            #         print(stderr, "$msg\n")
            #         error(msg)
            #     end
            # end
            # ccall(:jl_breakpoint, Cvoid, (Any,), slotptr)
            slot = UnsafeAtomics.load(slotptr)::MPCRQSlot{UInt32}
            (; index, safe, storage) = slot
            if !iszero(storage)
                if index == p_idx && slot.polarity !=′ polarity
                    threadindex = storage % Int
                    if x isa Waiter{T}
                        if inlinedata(crq)
                            y = _extract(T, storage)
                        else
                            y = crq.data[threadindex][itemindex]
                        end
                    else
                        y = crq.waiters[threadindex][itemindex]
                    end

                    newslot = MPCRQSlot{UInt32}(; safe, index = p_idx + R, polarity)
                    old = UnsafeAtomics.cas!(slotptr, slot, newslot)
                    if old == slot
                        if x isa Waiter{T}
                            return Some{T}(y)
                        else
                            return y
                        end
                    end
                else
                    newslot = MPCRQSlot(; safe = false, index, slot.polarity, storage)
                    old = UnsafeAtomics.cas!(slotptr, slot, newslot)
                    if old == slot
                        break
                    end
                end
                # If no `break` or `return`, reload the slot.
            else
                # Getting here means slot is empty.  Try to enqueue `x`.
                if x isa Waiter{T}
                    q = @atomic crq.data_idx
                else
                    q = @atomic crq.antidata_idx
                end
                q_idx = extract_index(q)
                if safe  # || q_idx ≤ p_idx
                    # if x isa Waiter{T}
                    #     GC.safepoint()
                    #     continue  # spin
                    # end
                    threadindex = Threads.threadid()
                    storage = threadindex % UInt32
                    if x isa Waiter{T}
                        crq.waiters[threadindex][itemindex] = x
                    else
                        if inlinedata(crq)
                            storage = _embed(UInt32, x)
                        else
                            crq.data[threadindex][itemindex] = x
                        end
                    end

                    @assert !iszero(storage)
                    newslot = MPCRQSlot(; safe = true, index = p_idx, polarity, storage)
                    old = UnsafeAtomics.cas!(slotptr, slot, newslot)
                    if old == slot
                        # if slot.index > p_idx
                        #     global CRQ = crq
                        # end
                        return MPCRQ_ENQUEUED
                        # else: reload the slot
                    end
                else
                    break
                end
            end
            GC.safepoint()  # cannot use yield
        end  # of the CAS loop on the slot
        # ccall(:jl_breakpoint, Cvoid, (Any,), nothing)

        # Reaching here means that a slot is successfully marked as unsafe.
        # Check if the CRQ is full and detect starvation.

        # Load the opposite index:
        if x isa Waiter{T}
            q = @atomic crq.data_idx
        else
            q = @atomic crq.antidata_idx
        end
        q_idx = extract_index(q)

        starvation_counter += 1
        if (p_idx - q_idx ≥ R) || starvation_counter > CRQ_STARVATION
            if x isa Waiter{T}
                set_closing!(crq, Val(:antidata_idx))
            else
                set_closing!(crq, Val(:data_idx))
            end
            if discovered_closing(crq) ≤ p_idx
                return MPCRQ_CLOSED
            end
        end
        GC.safepoint()  # cannot use yield
    end
end

# TODO: use bit-wise or?
function set_closing!(
    crq::IndirectMultiPolarityConcurrentRingQueueNode,
    ::Val{field},
) where {field}
    idx = getfield(crq, field, :sequentially_consistent)
    idx < 0 && return idx
    while true
        new = set_closing(idx)
        idx, success = replacefield!(
            crq,
            field,
            idx,
            new,
            :sequentially_consistent,
            :sequentially_consistent,
        )
        success && return new
        idx < 0 && return idx
    end
end

@noinline function discovered_closing(crq::IndirectMultiPolarityConcurrentRingQueueNode)
    closed_info = @atomic crq.closed_info
    closed_info_idx = extract_index(closed_info)
    isclosed = closed_info < 0
    isclosed && return closed_info_idx

    antidata_idx = set_closing!(crq, Val(:antidata_idx))
    data_idx = set_closing!(crq, Val(:data_idx))

    # @assert crq.antidata_idx < 0
    # @assert crq.data_idx < 0

    new = set_closing(max(extract_index(data_idx), extract_index(antidata_idx)))
    # @assert new < 0
    closed_info, success = @atomicreplace crq.closed_info closed_info => new
    # @assert Implementations.isclosed(crq)
    if success
        return extract_index(new)
    else
        return extract_index(closed_info)
    end
end

mutable struct DualLinkedConcurrentRingQueue{
    T,
    CRQ<:IndirectMultiPolarityConcurrentRingQueueNode{T},
}
    # _tag_pad::PadAfter32
    @atomic data::CRQ
    _data_pad::PadAfter64
    @atomic antidata::CRQ
    _antidata_pad::PadAfter64
    crqcache::ThreadLocalCache{CRQ}
    waitercache::ThreadLocalCache{Waiter{T}}
end

DualLinkedConcurrentRingQueue(; kwargs...) = DualLinkedConcurrentRingQueue{Any}(; kwargs...)
function DualLinkedConcurrentRingQueue{T}(; log2ringsize = 11) where {T}
    node = IndirectMultiPolarityConcurrentRingQueueNode{T}(log2ringsize)
    return DualLinkedConcurrentRingQueue(
        # _PadAfter32(),
        node,
        _PadAfter64(),
        node,
        _PadAfter64(),
        ThreadLocalCache{typeof(node)}(),
        ThreadLocalCache{Waiter{T}}(),
    )::DualLinkedConcurrentRingQueue{T}
end

Base.eltype(::Type{<:DualLinkedConcurrentRingQueue{T}}) where {T} = T

function Base.push!(lcrq::DualLinkedConcurrentRingQueue{T}, x) where {T}
    x = convert(eltype(lcrq), x)
    while true
        y = denqueue!(lcrq, x)
        if y isa Waiter{T}
            tryput!(y, x) || continue
        else
            @assert y === MPCRQ_ENQUEUED
        end
        return lcrq
    end
end

function Base.popfirst!(lcrq::DualLinkedConcurrentRingQueue{T}) where {T}
    w = let cached = _maybepop!(lcrq.waitercache)
        if cached === nothing
            Waiter{eltype(lcrq)}()
        else
            reinit!(something(cached))
        end
    end
    y = denqueue!(lcrq, w)
    if y === MPCRQ_ENQUEUED
        x = fetch(w)
    else
        trash!(w)
        x = something(y::Some{T})
    end
    trypush!(lcrq.waitercache, w)
    return x::T
end

# const CLOSED_X = Vector{Int}[Int[] for _ in 1:Threads.nthreads()]

"""
    denqueue!(lcrq::DualLCRQ{T}, x::T) -> MPCRQ_ENQUEUED or y::Waiter
    denqueue!(lcrq::DualLCRQ{T}, x::Waiter{T}) -> MPCRQ_ENQUEUED or y::Some{T}
    where DualLCRQ = DualLinkedConcurrentRingQueue
"""
function denqueue!(lcrq::DualLinkedConcurrentRingQueue{T}, x::Union{T,Waiter{T}}) where {T}
    while true
        crq = if x isa Waiter{T}
            @atomic lcrq.antidata
        else
            @atomic lcrq.data
        end
        y = denqueue!(crq, x)
        y === MPCRQ_CLOSED || return y

        # crq is closed (but may NOT be empty)
        next = @atomic crq.next
        if next !== nothing
            if x isa Waiter{T}
                @atomicreplace lcrq.antidata crq => next
            else
                @atomicreplace lcrq.data crq => next
            end
            # success && isempty(crq) && cache_crq!(lcrq, crq)
        else
            newcrq = make_newcrq!(lcrq, crq)
            # newcrq = similar(crq)
            y = denqueue!(newcrq, x)
            @assert y === MPCRQ_ENQUEUED
            _, success = @atomicreplace crq.next nothing => newcrq
            if success
                if x isa Waiter{T}
                    @atomicreplace lcrq.antidata crq => newcrq
                    return y
                else
                    @atomicreplace lcrq.data crq => newcrq
                    return y
                end
            else
                cache_crq!(lcrq, newcrq)
            end
        end
        GC.safepoint()  # yield?
    end
end

function make_newcrq!(lcrq::DualLinkedConcurrentRingQueue, crq)
    oldcrq = _maybepop!(lcrq.crqcache)
    if oldcrq === nothing
        return similar(crq)
    else
        return oldcrq
    end
end

function cache_crq!(lcrq::DualLinkedConcurrentRingQueue, crq)
    trypush!(lcrq.crqcache, empty!(crq))
    return
end

function nopop_foreach(f, crq::IndirectMultiPolarityConcurrentRingQueueNode)
    i = extract_index(@atomic crq.antidata_idx)
    while true
        slot = crq.ring[_mod1(i, crq)]
        (; index, polarity, storage) = slot
        if index == i && polarity ==′ DATA
            (; threadindex, itemindex) = storage
            x = crq.data[threadindex][itemindex]
            f(x)
            i += one(i)
        else
            break
        end
    end
end

function nopop_foreach(f, lcrq::DualLinkedConcurrentRingQueue)
    crq = @atomic lcrq.antidata
    while true
        nopop_foreach(f, crq)
        crq = @atomic crq.next
        crq === nothing && return
    end
end

datalength(crq::IndirectMultiPolarityConcurrentRingQueueNode) =
    extract_index(@atomic crq.data_idx) - extract_index(@atomic crq.antidata_idx)

isclosed(crq::IndirectMultiPolarityConcurrentRingQueueNode) = (@atomic crq.closed_info) < 0

function Base.show(io::IO, ::MIME"text/plain", lcrq::DualLinkedConcurrentRingQueue)
    data = @atomic lcrq.data
    antidata = @atomic lcrq.antidata
    if (@atomic antidata.next) === nothing && datalength(antidata) < 0
        node = data  # head of antidata
        polarity = ANTIDATA
    else
        node = antidata  # head of data
        polarity = DATA
    end
    n = 0
    ncrqs = 0
    while true
        n += datalength(node)
        ncrqs += 1
        node = @atomic node.next
        node === nothing && break
    end
    if polarity == DATA
        items = "$n data item" * (n > 1 ? "s" : "")
    else
        items = "$(-n) waiter" * (n < -1 ? "s" : "")
    end
    crqs = "$ncrqs node" * (ncrqs > 1 ? "s" : "")
    print(io, "Dual LCRQ: $items in $crqs")
end

function Base.show(
    io::IO,
    ::MIME"text/plain",
    crq::IndirectMultiPolarityConcurrentRingQueueNode,
)
    data_idx = extract_index(@atomic crq.data_idx)
    antidata_idx = extract_index(@atomic crq.antidata_idx)
    n = data_idx - antidata_idx
    if n == 0
        items = "empty"
    elseif n > 0
        items = "$n data item" * (n > 1 ? "s" : "")
    else
        items = "$(-n) waiter" * (n < -1 ? "s" : "")
    end
    status = isclosed(crq) ? "closed" : "open"
    print(
        io,
        "MP-CRQ: $items (status: $status, data_idx: $data_idx, antidata_idx: $antidata_idx)",
    )
end

function Base.NamedTuple(slot::MPCRQSlot)
    (; index, safe, polarity, storage) = slot
    return (; index, safe, polarity, storage)
end

function Base.show(io::IO, ::MIME"text/plain", slot::MPCRQSlot)
    if get(io, :typeinfo, Any) !== typeof(slot)
        show(io, MIME"text/plain"(), typeof(slot))
    end

    (; index, safe, polarity, storage) = slot
    nt = (;
        index,
        safe,
        polarity = Text(string(polarity)),
        storage = Text(
            sprint(
                show,
                MIME"text/plain"(),
                storage;
                context = IOContext(io, :typeinfo => typeof(storage)),
            ),
        ),
    )
    show(io, MIME"text/plain"(), nt)
end


text:
# Blender
export blender

#=
sel = utils.selection_get()
                            bpy.ops.view3d.select(location=(event.mouse_region_x, event.mouse_region_y))
                            sel1 = utils.selection_get()
                            if sel[0] != sel1[0] and sel1[0].type != 'MESH':
                                object = sel1[0]
                                target_slot = sel1[0].active_material_index
                                ui_props.has_hit = True
                            utils.selection_set(sel)

def mouse_raycast(context, mx, my):
    r = context.region
    rv3d = context.region_data
    coord = mx, my

    # get the ray from the viewport and mouse
    view_vector = view3d_utils.region_2d_to_vector_3d(r, rv3d, coord)
    ray_origin = view3d_utils.region_2d_to_origin_3d(r, rv3d, coord)
    ray_target = ray_origin + (view_vector * 1000000000)

    vec = ray_target - ray_origin

    has_hit, snapped_location, snapped_normal, face_index, object, matrix = bpy.context.scene.ray_cast(
        bpy.context.view_layer, ray_origin, vec)

    # rote = mathutils.Euler((0, 0, math.pi))
    randoffset = math.pi
    if has_hit:
        snapped_rotation = snapped_normal.to_track_quat('Z', 'Y').to_euler()
        up = Vector((0, 0, 1))
        props = bpy.context.scene.blenderkit_models
        if props.randomize_rotation and snapped_normal.angle(up) < math.radians(10.0):
            randoffset = props.offset_rotation_amount + math.pi + (
                    random.random() - 0.5) * props.randomize_rotation_amount
        else:
            randoffset = props.offset_rotation_amount  # we don't rotate this way on walls and ceilings. + math.pi
        # snapped_rotation.z += math.pi + (random.random() - 0.5) * .2

    else:
        snapped_rotation = mathutils.Quaternion((0, 0, 0, 0)).to_euler()

    snapped_rotation.rotate_axis('Z', randoffset)

    return has_hit, snapped_location, snapped_normal, snapped_rotation, face_index, object, matrix

=#

# BLR is a subtype of Python
parse_signature(::Val{:BLR}, sig::T) where {T} = parse_signature(Val(:PY), sig)
encode(::Val{:BLR}, t::Val{T}, c::IO, v) where {T} = encode(Val(:PY), t, c, v)
decode(::Val{:BLR}, t::Val{T}, c::IO) where {T} = decode(Val(:PY), t, c)
encode(ns::Val{:BLR}, t::Tuple{T1,T2,T3}, c::IO, v) where {T1,T2,T3} =
  begin
    encode(ns, T1(), c, v[1])
    encode(ns, T2(), c, v[2])
    encode(ns, T3(), c, v[3])
  end
decode(ns::Val{:BLR}, t::Tuple{T1,T2,T3}, c::IO) where {T1,T2,T3} =
  (decode(ns, T1(), c),
   decode(ns, T2(), c),
   decode(ns, T3(), c))
#=
# We need some additional Encoders
@encode_decode_as(:BLR, Val{:Entity}, Val{:size})
@encode_decode_as(:BLR, Val{:ObjectId}, Val{:size})
@encode_decode_as(:BLR, Val{:BIMLevel}, Val{:size})
@encode_decode_as(:BLR, Val{:FloorFamily}, Val{:size})
=#
@encode_decode_as(:BLR, Val{:Id}, Val{:size})
@encode_decode_as(:BLR, Val{:MatId}, Val{:size})

encode(::Val{:BLR}, t::Union{Val{:Point3d},Val{:Vector3d}}, c::IO, p) =
  encode(Val(:PY), Val(:float3), c, raw_point(p))
decode(::Val{:BLR}, t::Val{:Point3d}, c::IO) =
  xyz(decode(Val(:PY), Val(:float3), c)..., world_cs)
decode(::Val{:BLR}, t::Val{:Vector3d}, c::IO) =
  vxyz(decode(Val(:PY), Val(:float3), c)..., world_cs)

encode(ns::Val{:BLR}, t::Val{:Frame3d}, c::IO, v) = begin
  encode(ns, Val(:Point3d), c, v)
  t = v.cs.transform
  encode(Val(:PY), Val(:float3), c, (t[1,1], t[2,1], t[3,1]))
  encode(Val(:PY), Val(:float3), c, (t[1,2], t[2,2], t[3,2]))
  encode(Val(:PY), Val(:float3), c, (t[1,3], t[2,3], t[3,3]))
end

decode(ns::Val{:BLR}, t::Val{:Frame3d}, c::IO) =
  u0(cs_from_o_vx_vy_vz(
      decode(ns, Val(:Point3d), c),
      decode(ns, Val(:Vector3d), c),
      decode(ns, Val(:Vector3d), c),
      decode(ns, Val(:Vector3d), c)))

blender_api = @remote_functions :BLR """
def find_or_create_collection(name:str, active:bool, color:RGBA)->str:
def get_current_collection()->str:
def set_current_collection(name:str)->None:
def delete_all_shapes_in_collection(name:str)->None:
def delete_all_shapes()->None:
def delete_shape(name:Id)->None:
def select_shape(name:Id)->None:
def deselect_shape(name:Id)->None:
def deselect_all_shapes()->None:
def get_material(name:str)->MatId:
def get_blenderkit_material(ref:str)->MatId:
def new_material(name:str, diffuse_color:RGBA, specularity:float, roughness:float)->MatId:
def line(ps:List[Point3d], closed:bool, mat:MatId)->Id:
def spline(ps:List[Point3d], closed:bool, mat:MatId)->Id:
def nurbs(order:int, ps:List[Point3d], closed:bool, mat:MatId)->Id:
def mesh(verts:List[Point3d], edges:List[Tuple[int,int]], faces:List[List[int]], mat:MatId)->Id:
def objmesh(verts:List[Point3d], edges:List[Tuple[int,int]], faces:List[List[int]], smooth:bool, mat:MatId)->Id:
def trig(p1:Point3d, p2:Point3d, p3:Point3d, mat:MatId)->Id:
def quad(p1:Point3d, p2:Point3d, p3:Point3d, p4:Point3d, mat:MatId)->Id:
def quad_strip(ps:List[Point3d], qs:List[Point3d], smooth:bool, mat:MatId)->Id:
def quad_strip_closed(ps:List[Point3d], qs:List[Point3d], smooth:bool, mat:MatId)->Id:
def ngon(ps:List[Point3d], pivot:Point3d, smooth:bool, mat:MatId)->Id:
def polygon(ps:List[Point3d], mat:MatId)->Id:
def polygon_with_holes(pss:List[List[Point3d]], mat:MatId)->Id:
def quad_surface(ps:List[Point3d], nu:int, nv:int, closed_u:bool, closed_v:bool, smooth:bool, mat:MatId)->Id:
def circle(c:Point3d, v:Vector3d, r:float, mat:MatId)->Id:
def cuboid(verts:List[Point3d], mat:MatId)->Id:
def pyramid_frustum(bs:List[Point3d], ts:List[Point3d], smooth:bool, bmat:MatId, tmat:MatId, smat:MatId)->Id:
def sphere(center:Point3d, radius:float, mat:MatId)->Id:
def cone_frustum(b:Point3d, br:float, t:Point3d, tr:float, bmat:MatId, tmat:MatId, smat:MatId)->Id:
def box(p:Point3d, vx:Vector3d, vy:Vector3d, dx:float, dy:float, dz:float, mat:MatId)->Id:
def text(txt:str, p:Point3d, vx:Vector3d, vy:Vector3d, size:float)->Id:
def area_light(p:Point3d, v:Vector3d, size:float, color:RGBA, strength:float)->Id:
def sun_light(p:Point3d, v:Vector3d)->Id:
def light(p:Point3d, type:str)->Id:
def camera_from_view()->None:
def set_view(camera:Point3d, target:Point3d, lens:float)->None:
def get_view()->Tuple[Point3d, Point3d, float]:
def set_camera_view(camera:Point3d, target:Point3d, lens:float)->None:
def set_render_size(width:int, height:int)->None:
def set_render_path(filepath:str)->None:
def default_renderer()->None:
def cycles_renderer(samples:int, denoising:bool, motion_blur:bool, transparent:bool)->None:
def set_sun(latitude:float, longitude:float, elevation:float, year:int, month:int, day:int, time:float, UTC_zone:float, use_daylight_savings:bool)->None:
def set_sky(turbidity:float)->None:
"""

abstract type BLRKey end
const BLRId = Union{Int32,String} # Although shapes and materials are ints, layers are strings
const BLRIds = Vector{BLRId}
const BLRRef = NativeRef{BLRKey, BLRId}
const BLRRefs = Vector{BLRRef}
const BLREmptyRef = EmptyRef{BLRKey, BLRId}
const BLRUniversalRef = UniversalRef{BLRKey, BLRId}
const BLRUnionRef = UnionRef{BLRKey, BLRId}
const BLRSubtractionRef = SubtractionRef{BLRKey, BLRId}
const BLR = SocketBackend{BLRKey, BLRId}

const KhepriServerPath = Parameter(abspath(@__DIR__, "KhepriServer.py"))
export headless_blender
const headless_blender = Parameter(false)
const starting_blender = Parameter(false)

start_blender() =
  starting_blender() ?
    sleep(1) : # Just wait a little longer
    let blender_cmd = Sys.iswindows() ?
      	  joinpath(readdir("C:/Program Files/Blender Foundation/", join=true)[1], "blender.exe") :
    	  "blender"
	  starting_blender(true)
      run(detach(headless_blender() ?
            `$(blender_cmd) -noaudio --background --python $(KhepriServerPath())` :
      	    `$(blender_cmd) --python $(KhepriServerPath())`),
    	  wait=false)
    end

KhepriBase.retry_connecting(b::BLR) =
  (@info("Starting $(b.name)."); start_blender(); sleep(2))

KhepriBase.after_connecting(b::BLR) =
  begin
	starting_blender(false)
	#set_material(blender, material_basic, )
	set_material(blender, material_metal, "asset_base_id:f1774cb0-b679-46b4-879e-e7223e2b4b5f asset_type:material")
	#set_material(blender, material_glass, "asset_base_id:ee2c0812-17f5-40d4-992c-68c5a66261d7 asset_type:material")
	set_material(blender, material_glass, "asset_base_id:ffa3c281-6184-49d8-b05e-8c6e9fe93e68 asset_type:material")
	set_material(blender, material_wood, "asset_base_id:d5097824-d5a1-4b45-ab5b-7b16bdc5a627 asset_type:material")
	#set_material(blender, material_concrete, "asset_base_id:0662b3bf-a762-435d-9407-e723afd5eafc asset_type:material")
	set_material(blender, material_concrete, "asset_base_id:df1161da-050c-4638-b376-38ced992ec18 asset_type:material")
	set_material(blender, material_plaster, "asset_base_id:c674137d-cfae-45f1-824f-e85dc214a3af asset_type:material")

	#set_material(blender, material_grass, "asset_base_id:97b171b4-2085-4c25-8793-2bfe65650266 asset_type:material")
	#set_material(blender, material_grass, "asset_base_id:7b05be22-6bed-4584-a063-d0e616ddea6a asset_type:material")
	set_material(blender, material_grass, "asset_base_id:b4be2338-d838-433b-9f0d-2aa9b97a0a8a asset_type:material")
  end

const blender = BLR("Blender", blender_port, blender_api)

KhepriBase.has_boolean_ops(::Type{BLR}) = HasBooleanOps{false}()

KhepriBase.backend(::BLRRef) = blender
KhepriBase.void_ref(b::BLR) = BLRRef(-1 % Int32)

# Primitives

KhepriBase.b_line(b::BLR, ps, mat) =
  @remote(b, line(ps, false, mat))

KhepriBase.b_polygon(b::BLR, ps, mat) =
	@remote(b, line(ps, true, mat))

KhepriBase.b_nurbs_curve(b::BLR, order, ps, knots, weights, closed, mat) =
  @remote(b, nurbs(order, ps, closed, mat))

KhepriBase.b_trig(b::BLR, p1, p2, p3, mat) =
  @remote(b, trig(p1, p2, p3, mat))

KhepriBase.b_quad(b::BLR, p1, p2, p3, p4, mat) =
	@remote(b, quad(p1, p2, p3, p4, mat))

KhepriBase.b_ngon(b::BLR, ps, pivot, smooth, mat) =
	@remote(b, ngon(ps, pivot, smooth, mat))

KhepriBase.b_quad_strip(b::BLR, ps, qs, smooth, mat) =
  @remote(b, quad_strip(ps, qs, smooth, mat))

KhepriBase.b_quad_strip_closed(b::BLR, ps, qs, smooth, mat) =
  @remote(b, quad_strip_closed(ps, qs, smooth, mat))

KhepriBase.b_surface_polygon(b::BLR, ps, mat) =
  @remote(b, polygon(ps, mat))

KhepriBase.b_surface_polygon_with_holes(b::BLR, ps, qss, mat) =
  @remote(b, polygon_with_holes([ps, qss...], mat))

KhepriBase.b_surface_circle(b::BLR, c, r, mat) =
  @remote(b, circle(c, vz(1, c.cs), r, mat))

KhepriBase.b_surface_grid(b::BLR, ptss, closed_u, closed_v, smooth_u, smooth_v, mat) =
  let (nu, nv) = size(ptss)
	smooth_u && smooth_v ?
	  @remote(b, quad_surface(vcat(ptss...), nu, nv, closed_u, closed_v, true, mat)) :
	  smooth_u ?
	  	(closed_u ?
          vcat([b_quad_strip_closed(b, ptss[:,i], ptss[:,i+1], true, mat) for i in 1:nv-1],
	           closed_v ? [b_quad_strip_closed(b, ptss[:,end], ptss[:,1], true, mat)] : []) :
	      vcat([b_quad_strip(b, ptss[:,i], ptss[:,i+1], true, mat) for i in 1:nv-1],
	           closed_v ? [b_quad_strip(b, ptss[:,end], ptss[:,1], true, mat)] : [])) :
 	    (closed_v ?
           vcat([b_quad_strip_closed(b, ptss[i,:], ptss[i+1,:], smooth_v, mat) for i in 1:nu-1],
  	         	closed_u ? [b_quad_strip_closed(b, ptss[end,:], ptss[1,:], smooth_v, mat)] : []) :
  	       vcat([b_quad_strip(b, ptss[i,:], ptss[i+1,:], smooth_v, mat) for i in 1:nu-1],
  	          	closed_u ? [b_quad_strip(b, ptss[end,:], ptss[1,:], smooth_v, mat)] : []))
  end

KhepriBase.b_generic_pyramid_frustum(b::BLR, bs, ts, smooth, bmat, tmat, smat) =
  @remote(b, pyramid_frustum(bs, ts, smooth, bmat, tmat, smat))

KhepriBase.b_cone(b::BLR, cb, r, h, bmat, smat) =
  @remote(b, cone_frustum(cb, r, add_z(cb, h), 0, bmat, bmat, smat))

KhepriBase.b_cone_frustum(b::BLR, cb, rb, h, rt, bmat, tmat, smat) =
  @remote(b, cone_frustum(cb, rb, add_z(cb, h), rt, bmat, tmat, smat))

KhepriBase.b_cylinder(b::BLR, cb, r, h, bmat, tmat, smat) =
  @remote(b, cone_frustum(cb, r, add_z(cb, h), r, bmat, tmat, smat))

KhepriBase.b_cuboid(b::BLR, pb0, pb1, pb2, pb3, pt0, pt1, pt2, pt3, mat) =
  @remote(b, cuboid([pb0, pb1, pb2, pb3, pt0, pt1, pt2, pt3], mat))

KhepriBase.b_sphere(b::BLR, c, r, mat) =
  @remote(b, sphere(c, r, mat))

# Materials

KhepriBase.b_get_material(b::BLR, ref) =
  get_blender_material(b, ref)

#=
Important source of materials:

1. Activate Blender's blenderkit addon:
https://www.blenderkit.com/get-blenderkit/

2. Browse BlenderKit's material database:
https://www.blenderkit.com/asset-gallery?query=category_subtree:material

3. Select material and copy reference, e.g.:
asset_base_id:ced25dc0-d461-42f7-aa03-85cb88f671a1 asset_type:material

4. Install material and retrive its id with:
b_get_material(blender, "asset_base_id:ced25dc0-d461-42f7-aa03-85cb88f671a1 asset_type:material")
=#

get_blender_material(b, ref::Nothing) =
  void_ref(b)

get_blender_material(b, ref::AbstractString) =
  startswith(ref, "asset_base_id") ?
    @remote(b, get_blenderkit_material(ref)) :
    @remote(b, get_material(ref))

KhepriBase.b_new_material(b::BLR, path, color, specularity, roughness, transmissivity, transmitted_specular) =
  @remote(b, new_material(path, convert(RGBA, color), specularity, roughness))

# BIM

#=

Default families

=#
export blender_family_materials
blender_family_materials(m1, m2=m1, m3=m2, m4=m3) = (materials=(m1, m2, m3, m4), )

KhepriBase.b_layer(b::BLR, name, active, color) =
  @remote(b, find_or_create_collection(name, active, color))
KhepriBase.b_current_layer(b::BLR) =
  @remote(b, get_current_collection())
KhepriBase.b_current_layer(b::BLR, layer) =
  @remote(b, set_current_collection(layer))
KhepriBase.b_all_shapes_in_layer(b::BLR, layer) =
  @remote(b, all_shapes_in_collection(layer))
KhepriBase.b_delete_all_shapes_in_layer(b::BLR, layer) =
  @remote(b, delete_all_shapes_in_collection(layer))

realize(b::BLR, s::EmptyShape) =
  BLREmptyRef()
realize(b::BLR, s::UniversalShape) =
  BLRUniversalRef()

KhepriBase.b_set_view(b::BLR, camera::Loc, target::Loc, lens::Real, aperture::Real) =
  begin
  	@remote(b, set_view(camera, target, lens))
  	@remote(b, set_camera_view(camera, target, lens))
  end

KhepriBase.b_get_view(b::BLR) =
  @remote(b, get_view())

KhepriBase.b_zoom_extents(b::BLR) = @remote(b, ZoomExtents())

KhepriBase.b_set_view_top(b::BLR) = @remote(b, ViewTop())

KhepriBase.b_realistic_sky(b::BLR, date, latitude, longitude, elevation, meridian, turbidity, withsun) =
  begin
	@remote(b, set_sun(latitude, longitude, elevation, year(date), month(date), day(date), hour(date)+minute(date)/60, meridian, false))
	@remote(b, set_sky(turbidity)) #Add withsun
  end

KhepriBase.b_delete_ref(b::BLR, r::BLRId) =
  @remote(b, delete_shape(r))

KhepriBase.b_delete_all_refs(b::BLR) =
  @remote(b, delete_all_shapes())

#=
backend_set_length_unit(b::BLR, unit::String) = @remote(b, SetLengthUnit(unit))

# Dimensions

const BLRDimensionStyles = Dict(:architectural => "_ARCHTICK", :mechanical => "")

backend_dimension(b::BLR, p0::Loc, p1::Loc, p::Loc, scale::Real, style::Symbol) =
    @remote(b, CreateAlignedDimension(p0, p1, p,
        scale,
        BLRDimensionStyles[style]))

backend_dimension(b::BLR, p0::Loc, p1::Loc, sep::Real, scale::Real, style::Symbol) =
    let v = p1 - p0
        angle = pol_phi(v)
        dimension(p0, p1, add_pol(p0, sep, angle + pi/2), scale, style, b)
    end

# Blocks

realize(b::BLR, s::Block) =
    @remote(b, CreateBlockFromShapes(s.name, collect_ref(s.shapes)))

realize(b::BLR, s::BlockInstance) =
    @remote(b, CreateBlockInstance(
        collect_ref(s.block)[1],
        center_scaled_cs(s.loc, s.scale, s.scale, s.scale)))
=#
#=

# Manual process
@time for i in 1:1000 for r in 1:10 circle(x(i*10), r) end end

# Create block...
Khepri.create_block("Foo", [circle(radius=r) for r in 1:10])

# ...and instantiate it
@time for i in 1:1000 Khepri.instantiate_block("Foo", x(i*10), 0) end

=#
#=
# Lights
KhepriBase.b_pointlight(b::BLR, loc::Loc, color::RGB, range::Real, intensity::Real) =
  # HACK: Fix this
  @remote(b, SpotLight(loc, intensity, range, loc+vz(-1)))

KhepriBase.b_spotlight(b::BLR, loc::Loc, dir::Vec, hotspot::Real, falloff::Real) =
    @remote(b, SpotLight(loc, hotspot, falloff, loc + dir))

KhepriBase.b_ieslight(b::BLR, file::String, loc::Loc, dir::Vec, alpha::Real, beta::Real, gamma::Real) =
    @remote(b, IESLight(file, loc, loc + dir, vxyz(alpha, beta, gamma)))

# User Selection

backend_shape_from_ref(b::BLR, r) =
  let c = connection(b),
      code = @remote(b, ShapeCode(r)),
      ref = DynRefs(b=>BLRRef(r))
    if code == 1 # Point
        point(@remote(b, PointPosition(r)),
              ref=ref)
    elseif code == 2
        circle(loc_from_o_vz(@remote(b, CircleCenter(r)), @remote(b, CircleNormal(r))),
               @remote(b, CircleRadius(r)),
               ref=ref)
    elseif 3 <= code <= 6
        line(@remote(b, LineVertices(r)),
             ref=ref)
    elseif code == 7
        let tans = @remote(b, SplineTangents(r))
            if length(tans[1]) < 1e-20 && length(tans[2]) < 1e-20
                closed_spline(@remote(b, SplineInterpPoints(r))[1:end-1],
                              ref=ref)
            else
                spline(@remote(b, SplineInterpPoints(r)), tans[1], tans[2],
                       ref=ref)
            end
        end
    elseif code == 9
        let start_angle = mod(@remote(b, ArcStartAngle(r)), 2pi),
            end_angle = mod(@remote(b, ArcEndAngle(r)), 2pi)
            arc(loc_from_o_vz(@remote(b, ArcCenter(r)), @remote(b, ArcNormal(r))),
                @remote(b, ArcRadius(r)), start_angle, mod(end_angle - start_angle, 2pi),
                ref=ref)
        #=    if end_angle > start_angle
                arc(maybe_loc_from_o_vz(@remote(b, ArcCenter(r)), @remote(b, ArcNormal(r))),
                    @remote(b, ArcRadius(r)), start_angle, end_angle - start_angle,
                    ref=ref)
            else
                arc(maybe_loc_from_o_vz(@remote(b, ArcCenter(r)), @remote(b, ArcNormal(r))),
                    @remote(b, ArcRadius(r)), end_angle, start_angle - end_angle,
                    ref=ref)
            end=#
        end
    elseif code == 10
        let str = @remote(b, TextString(r)),
            height = @remote(b, TextHeight(r)),
            loc = @remote(b, TextPosition(r))
            text(str, loc, height, ref=ref)
        end
    elseif code == 11
        let str = @remote(b, MTextString(r)),
            height = @remote(b, MTextHeight(r)),
            loc = @remote(b, MTextPosition(r))
            text(str, loc, height, ref=ref)
        end
    elseif code == 16
        let pts = @remote(b, MeshVertices(r)),
            (type, n, m, n_closed, m_closed) = @remote(b, PolygonMeshData(r))
            surface_grid(reshape(pts, (n, m)), n_closed == 1, m_closed == 1, ref=ref)
        end
    elseif 12 <= code <= 14
        surface(Shapes1D[], ref=ref)
    elseif 103 <= code <= 106
        polygon(@remote(b, LineVertices(r)),
                ref=ref)
    elseif code == 107
        closed_spline(@remote(b, SplineInterpPoints(r))[1:end-1],
                      ref=ref)
    else
        #unknown(ref=ref)
        unknown(r, ref=ref) # To force copy
        #error("Unknown shape with code $(code)")
    end
  end
#
=#
#=
In case we need to realize an Unknown shape, we just copy it
=#
#=
realize(b::BLR, s::Unknown) =
    @remote(b, Copy(s.baseref))

backend_select_position(b::BLR, prompt::String) =
  begin
    @info "$(prompt) on the $(b) backend."
    let ans = @remote(b, GetPosition(prompt))
      length(ans) > 0 ? ans[1] : nothing
    end
  end

backend_select_positions(b::BLR, prompt::String) =
  let sel() =
    let p = select_position(prompt, b)
      if p == nothing
        []
      else
        [p, sel()...]
      end
    end
    sel()
  end

# HACK: The next operations should receive a set of shapes to avoid re-creating already existing shapes

backend_select_point(b::BLR, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetPoint)
backend_select_points(b::BLR, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetPoints)

backend_select_curve(b::BLR, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetCurve)
backend_select_curves(b::BLR, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetCurves)

backend_select_surface(b::BLR, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetSurface)
backend_select_surfaces(b::BLR, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetSurfaces)

backend_select_solid(b::BLR, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetSolid)
backend_select_solids(b::BLR, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetSolids)

backend_select_shape(b::BLR, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetShape)
backend_select_shapes(b::BLR, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetShapes)

backend_captured_shape(b::BLR, handle) =
  backend_shape_from_ref(b, @remote(b, GetShapeFromHandle(handle)))
backend_captured_shapes(b::BLR, handles) =
  map(handles) do handle
      backend_shape_from_ref(b, @remote(b, GetShapeFromHandle(handle)))
  end

backend_generate_captured_shape(b::BLR, s::Shape) =
    println("captured_shape(blender, $(@remote(b, GetHandleFromShape(ref(s).value))))")
backend_generate_captured_shapes(b::BLR, ss::Shapes) =
  begin
    print("captured_shapes(blender, [")
    for s in ss
      print(@remote(b, GetHandleFromShape(ref(s).value)))
      print(", ")
    end
    println("])")
  end

# Register for notification

backend_register_shape_for_changes(b::BLR, s::Shape) =
    let conn = connection(b)
        @remote(b, RegisterForChanges(ref(s).value))
        @remote(b, DetectCancel())
        s
    end

backend_unregister_shape_for_changes(b::BLR, s::Shape) =
    let conn = connection(b)
        @remote(b, UnregisterForChanges(ref(s).value))
        @remote(b, UndetectCancel())
        s
    end

backend_waiting_for_changes(b::BLR, s::Shape) =
    ! @remote(b, WasCanceled())

backend_changed_shape(b::BLR, ss::Shapes) =
    let conn = connection(b)
        changed = []
        while length(changed) == 0 && ! @remote(b, WasCanceled())
            changed =  @remote(b, ChangedShape())
            sleep(0.1)
        end
        if length(changed) > 0
            backend_shape_from_ref(b, changed[1])
        else
            nothing
        end
    end

# HACK: This should be filtered on the plugin, not here.
backend_all_shapes(b::BLR) =
  Shape[backend_shape_from_ref(b, r)
        for r in filter(r -> @remote(b, ShapeCode(r)) != 0, @remote(b, GetAllShapes()))]

backend_all_shapes_in_layer(b::BLR, layer) =
  Shape[backend_shape_from_ref(b, r) for r in @remote(b, GetAllShapesInLayer(layer))]

=#

KhepriBase.b_highlight_ref(b::BLR, r::BLRId) =
  @remote(b, select_shape(r))

KhepriBase.b_unhighlight_ref(b::BLR, r::BLRId) =
  @remote(b, deselect_shape(r))

KhepriBase.b_unhighlight_all_refs(b::BLR) =
  @remote(b, deselect_all_shapes())
#=
backend_pre_selected_shapes_from_set(ss::Shapes) =
  length(ss) == 0 ? [] : pre_selected_shapes_from_set(ss, backend(ss[1]))

# HACK: This must be implemented for all backends
backend_pre_selected_shapes_from_set(ss::Shapes, b::Backend) = []

backend_pre_selected_shapes_from_set(b::BLR, ss::Shapes) =
  let refs = map(id -> @remote(b, GetHandleFromShape(id)), @remote(b, GetPreSelectedShapes()))
    filter(s -> @remote(b, GetHandleFromShape(ref(s).value)) in refs, ss)
  end
backend_disable_update(b::BLR) =
  @remote(b, DisableUpdate())

backend_enable_update(b::BLR) =
  @remote(b, EnableUpdate())

# Render

#render exposure: [-3, +3] -> [-6, 21]
convert_render_exposure(b::BLR, v::Real) = -4.05*v + 8.8
#render quality: [-1, +1] -> [+1, +50]
convert_render_quality(b::BLR, v::Real) = round(Int, 25.5 + 24.5*v)

=#

KhepriBase.b_render_view(b::BLR, path::String) =
  let (camera, target, lens) = @remote(b, get_view())
    @remote(b, set_camera_view(camera, target, lens))
    @remote(b, set_render_size(render_width(), render_height()))
    @remote(b, set_render_path(path))
  	@remote(b, cycles_renderer(200, true, false, false))
  end

#=

export blender_command
blender_command(s::String) =
  @remote(blender, Command("$(s)\n"))
=#


text:
using Documenter

push!(LOAD_PATH,  "../../src")
push!(LOAD_PATH,  "../../serializers")

using SearchLight, SearchLight.Callbacks, SearchLight.Configuration, SearchLight.Exceptions
using SearchLight.Generator.FileTemplates, SearchLight.Generator, SearchLight.Migrations, SearchLight.QueryBuilder
using SearchLight.Relationships, SearchLight.Serializer, SearchLight.Transactions, SearchLight.Validation
using SearchLight.Serializer.JsonSerializer

makedocs(
    sitename = "SearchLight - Concise, secure, cross-platform query builder and ORM for Julia",
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "Home" => "index.md",
        "SearchLight API" => [
          "Callbacks" => "api/callbacks.md",
          "Configuration" => "api/configuration.md",
          "Exceptions" => "api/exceptions.md",
          "FileTemplates" => "api/filetemplates.md",
          "Generator" => "api/generator.md",
          "Migrations" => "api/migrations.md",
          "ModelTypes" => "api/modeltypes.md",
          "QueryBuilder" => "api/querybuilder.md",
          "Relationships" => "api/relationships.md",
          "SearchLight" => "api/searchlight.md",
          "Serializer" => "api/serializer.md",
          "Serializers" => [
            "JsonSerializer" => "api/serializers/json.md",
          ],
          "Transactions" => "api/transactions.md",
          "Validation" => "api/validation.md",
        ]
    ],
)

deploydocs(
  repo = "github.com/GenieFramework/SearchLight.jl.git",
)


text:
# Latafat, Patrinos, "Asymmetric forward–backward–adjoint splitting for
# solving monotone inclusions involving three operators", Computational
# Optimization and Applications, vol. 68, no. 1, pp. 57-93 (2017).
#
# Latafat, Patrinos, "Primal-dual proximal algorithms for structured convex
# optimization: a unifying framework", In Large-Scale and Distributed
# Optimization, Giselsson and Rantzer, Eds. Springer International Publishing,
# pp. 97–120 ( 2018).
#
# Chambolle, Pock, "A First-Order Primal-Dual Algorithm for Convex Problems
# with Applications to Imaging", Journal of Mathematical Imaging and Vision,
# vol. 40, no. 1, pp. 120-145 (2011).
#
# Condat, "A primal–dual splitting method for convex optimization
# involving Lipschitzian, proximable and linear composite terms",
# Journal of Optimization Theory and Applications, vol. 158, no. 2,
# pp 460-479 (2013).
#
# Vũ, "A splitting algorithm for dual monotone inclusions involving
# cocoercive operators", Advances in Computational Mathematics, vol. 38, no. 3,
# pp. 667-681 (2013).

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

struct AFBA_iterable{R, Tx, Ty, Tf, Tg, Th, Tl, TL}
    f::Tf
    g::Tg
    h::Th
    l::Tl
    L::TL
    x0::Tx
    y0::Ty
    theta::R
    mu::R
    lambda::R
    gamma1::R
    gamma2::R
end

Base.IteratorSize(::Type{<:AFBA_iterable}) = Base.IsInfinite()

struct AFBA_state{Tx, Ty}
    x::Tx
    y::Ty
    xbar::Tx
    ybar::Ty
    gradf::Tx
    gradl::Ty
    FPR_x::Tx
    FPR_y::Ty
    temp_x::Tx
    temp_y::Ty
end

AFBA_state(iter::AFBA_iterable) = AFBA_state(
    copy(iter.x0), copy(iter.y0),
    zero(iter.x0), zero(iter.y0),
    zero(iter.x0), zero(iter.y0),
    zero(iter.x0), zero(iter.y0),
    zero(iter.x0), zero(iter.y0)
)

function Base.iterate(iter::AFBA_iterable, state::AFBA_state=AFBA_state(iter))
    # perform xbar-update step
    gradient!(state.gradf, iter.f, state.x)
    mul!(state.temp_x, iter.L', state.y)
    state.temp_x .+= state.gradf
    state.temp_x .*= -iter.gamma1
    state.temp_x .+= state.x
    prox!(state.xbar, iter.g, state.temp_x, iter.gamma1)

    # perform ybar-update step
    gradient!(state.gradl, Conjugate(iter.l), state.y)
    state.temp_x .= iter.theta .* state.xbar .+ (1-iter.theta) .* state.x
    mul!(state.temp_y, iter.L, state.temp_x)
    state.temp_y .-= state.gradl
    state.temp_y .*= iter.gamma2
    state.temp_y .+= state.y
    prox!(state.ybar, Conjugate(iter.h), state.temp_y, iter.gamma2)

    # the residues
    state.FPR_x .= state.xbar .- state.x
    state.FPR_y .= state.ybar .- state.y

    # perform x-update step
    state.temp_y .= (iter.mu * (2-iter.theta) * iter.gamma1) .* state.FPR_y
    mul!(state.temp_x, iter.L', state.temp_y)
    state.x .+= iter.lambda .* (state.FPR_x .- state.temp_x)

    # perform y-update step
    state.temp_x .= ((1-iter.mu) * (2-iter.theta) * iter.gamma2) .* state.FPR_x
    mul!(state.temp_y, iter.L, state.temp_x)
    state.y .+= iter.lambda .* (state.FPR_y .+ state.temp_y)

    return state, state
end

function AFBA_default_stepsizes(L, h, theta::R, mu::R, beta_f::R, beta_l::R) where {R <: Real}
    # lambda = 1
    if isa(h, Zero)
        gamma1 = R(1.99)/beta_f
        gamma2 = R(1) # does not matter
    else
        par = R(5) # scaling parameter for comparing Lipschitz constants and \|L\|
        par2 = R(100)   # scaling parameter for α
        alpha = R(1)
        nmL = R(opnorm(L))
        if theta == 2 # default stepsize for Vu-Condat
            if nmL > par * max(beta_l, beta_f)
                alpha = R(1)
            elseif beta_f > par*beta_l
                alpha = par2*nmL/beta_f
            elseif beta_l > par*beta_f
                alpha = beta_l/(par2*nmL)
            end
            gamma1 = R(1)/(beta_f/2+nmL/alpha)
            gamma2 = R(0.99)/(beta_l/2+nmL*alpha)
        elseif theta == 1 && mu == 1 # SPCA
            if nmL > par2*beta_l # for the case beta_f = 0
                alpha = R(1)
            elseif beta_l > par*beta_f
                alpha =  beta_l/(par2*nmL)
            end
            gamma1 = beta_f > 0 ?  R(1.99)/beta_f : R(1)/(nmL/alpha)
            gamma2 = R(0.99)/(beta_l/2 + gamma1*nmL^2)
        elseif theta == 0 && mu == 1 # PPCA
            temp = R(3)
            if beta_f == 0
                nmL *= sqrt(temp)
                if nmL > par *beta_l
                    alpha = R(1)
                else
                    alpha = beta_l/(par2*nmL)
                end
                gamma1 = R(1)/(beta_f/2+nmL/alpha)
                gamma2 = R(0.99)/(beta_l/2+nmL*alpha)
            else
                if nmL > par * max(beta_l, beta_f)
                    alpha = R(1)
                elseif beta_f > par*beta_l
                    alpha = par2*nmL/beta_f
                elseif beta_l > par*beta_f
                    alpha = beta_l/(par2*nmL)
                end
                xi = 1+ 2*nmL/(nmL+alpha*beta_f/2)
                gamma1 = R(1)/(beta_f/2+nmL/alpha)
                gamma2 = R(0.99)/(beta_l/2+xi*nmL*alpha)
            end
        elseif mu == 0 # SDCA & PDCA
            temp = theta^2-3*theta+3
            if beta_l == R(0)
                nmL *= sqrt(temp)
                 if nmL > par *beta_f
                    alpha = R(1)
                else
                    alpha = par2*nmL/beta_f
                end
                gamma1 = R(1)/(beta_f/2+nmL/alpha)
                gamma2 = R(0.99)/(beta_l/2+nmL*alpha)
            else
                if nmL > par * max(beta_l, beta_f)
                    alpha = R(1)
                elseif beta_f > par*beta_l
                    alpha = par2*nmL/beta_f
                elseif beta_l > par*beta_f
                    alpha = beta_l/(par2*nmL)
                end
                eta = 1+ (temp-1)*alpha*nmL/(alpha*nmL+beta_l/2)
                gamma1 = R(1)/(beta_f/2+eta*nmL/alpha)
                gamma2 = R(0.99)/(beta_l/2+nmL*alpha)
            end
        elseif theta == 0 && mu == 0.5 # PPDCA
            if beta_l == 0 || beta_f == 0
                 if nmL > par * max(beta_l, beta_f)
                    alpha = R(1)
                elseif beta_f > par*beta_l
                    alpha = par2*nmL/beta_f
                elseif beta_l > par*beta_f
                    alpha = beta_l/(par2*nmL)
                end
            else
                alpha = sqrt(beta_l/beta_f)/2
            end
            gamma1 = R(1)/(beta_f/2+nmL/alpha)
            gamma2 = R(0.99)/(beta_l/2+nmL*alpha)
        else
            error("this choice of theta and mu is not supported!")
        end
    end

    return gamma1, gamma2
end

# Solver

struct AFBA{R}
    gamma1::Maybe{R}
    gamma2::Maybe{R}
    theta::R
    mu::R
    lambda::R
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int

    function AFBA{R}(; gamma1::Maybe{R}=nothing, gamma2::Maybe{R}=nothing,
        theta::R=R(1), mu::R=R(1), lambda::R=R(1), maxit::Int=10000,
        tol::R=R(1e-5), verbose::Bool=false, freq::Int=100
    ) where R
        @assert gamma1 === nothing || gamma1 > 0
        @assert gamma2 === nothing || gamma2 > 0
        @assert theta >= 0
        @assert 0 <= mu <= 1
        @assert lambda > 0
        @assert maxit > 0
        @assert tol > 0
        @assert freq > 0
        new(gamma1, gamma2, theta, mu, lambda, maxit, tol, verbose, freq)
    end
end

function (solver::AFBA{R})(x0::AbstractArray{C}, y0::AbstractArray{C};
    f=Zero(), g=Zero(), h=Zero(), l=IndZero(), L=I, beta_f::R=R(0), beta_l::R=R(0)
) where {R, C <: Union{R, Complex{R}}}

    stop(state::AFBA_state) = norm(state.FPR_x, Inf) + norm(state.FPR_y, Inf) <= solver.tol
    disp((it, state)) = @printf(
        "%6d | %7.4e\n",
        it, norm(state.FPR_x, Inf) + norm(state.FPR_y, Inf)
    )

    if isa(h, Zero) # case h(Lx) equiv 0
        L = R(0)*I
        y0 = zero(x0)
    end

    lambda = solver.lambda
    if solver.gamma1 === nothing || solver.gamma2 === nothing
        if solver.lambda != 1
            @warn "default stepsizes are not supported with this choice of lamdba,"*
            " reverted to default lambda"
            lambda = R(1)
        end
        gamma1, gamma2 = AFBA_default_stepsizes(
            L, h, solver.theta, solver.mu, beta_f, beta_l
        )
    else
        gamma1, gamma2 = solver.gamma1, solver.gamma2
    end

    iter = AFBA_iterable(f, g, h, l, L, x0, y0,
        solver.theta, solver.mu, lambda, gamma1, gamma2
    )
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose iter = tee(sample(iter, solver.freq), disp) end

    num_iters, state_final = loop(iter)

    return state_final.xbar, state_final.ybar, num_iters
end

# Outer constructors

"""
    AFBA([gamma1, gamma2, theta, mu, lambda, maxit, tol, verbose, freq])

Instantiate the asymmetric forward-backward-adjoing algorithm (AFBA, see [1])
for solving convex optimization problems of the form

    minimize f(x) + g(x) + (h □ l)(L x),

where `f` is smooth, `g` and `h` are possibly nonsmooth and `l` is strongly
convex. Symbol `□` denotes the infimal convolution, and `L` is a linear mapping.
If `solver = AFBA(args...)`, then the above problem is solved with

    solver(x0, y0; [f, g, h, l, L, beta_f, beta_l])

Points `x0` and `y0` are the initial primal and dual iterates, respectively.
If unspecified, functions `f`, `g`, and `h` default to the identically zero
function, `l` defaults to the indicator of the set `{0}`, and `L` defaults to
the identity. Important keyword arguments, in case `f` and `l` are set, are:

* `beta_f`: Lipschitz constant of gradient of `f` (default: zero)
* `beta_l`: Lipschitz constant of gradient of the conjugate of `l` (default: zero)

These are used to determine the algorithm default stepsizes, `gamma1` and `gamma2`,
in case they are not directly specified.

Optional keyword arguments are:

* `gamma1`: stepsize corresponding to the primal updates (default: see [1] for each case)
* `gamma2`: stepsize corresponding to the dual updates (default: see [1] for each case)
* `theta`: nonnegative algorithm parameter (default: `1.0`)
* `mu`: algorithm parameter in the range [0,1] (default: `1.0`)
* `tol`: primal-dual termination tolerance (default: `1e-5`)
* `maxit`: maximum number of iterations (default: `10000`)
* `verbose`, verbosity level (default: `1`)
* `verbose_freq`, verbosity frequency for `verbose = 1` (default: `100`)

The iterator implements Algorithm 3 of [1] with constant stepsize (α_n=λ)
for several prominant special cases:
1) θ = 2          ==>   Corresponds to the Vu-Condat Algorithm [2,3].
2) θ = 1, μ=1
3) θ = 0, μ=1
4) θ ∈ [0,∞), μ=0

See [2, Section 5.2] and [1, Figure 1] for stepsize conditions, special cases,
and relation to other algorithms.

References:

[1] Latafat, Patrinos, "Asymmetric forward–backward–adjoint splitting for
solving monotone inclusions involving three operators", Computational
Optimization and Applications, vol. 68, no. 1, pp. 57-93 (2017).

[2] Latafat, Patrinos, "Primal-dual proximal algorithms for structured convex
optimization : a unifying framework", In Large-Scale and Distributed
Optimization, Giselsson and Rantzer, Eds. Springer International Publishing,
pp. 97–120 ( 2018).

[3] Condat, "A primal–dual splitting method for convex optimization
involving Lipschitzian, proximable and linear composite terms",
Journal of Optimization Theory and Applications, vol. 158, no. 2,
pp 460-479 (2013).

[4] Vũ, "A splitting algorithm for dual monotone inclusions involving
cocoercive operators", Advances in Computational Mathematics, vol. 38, no. 3,
pp. 667-681 (2013).
"""
AFBA(::Type{R}; kwargs...) where R = AFBA{R}(; kwargs...)
AFBA(; kwargs...) = AFBA(Float64; kwargs...)

"""
    VuCondat([gamma1, gamma2, theta, mu, lambda, maxit, tol, verbose, freq])

Instantiate the Vû-Condat splitting algorithm (see [2, 3])
for solving convex optimization problems of the form

    minimize f(x) + g(x) + (h □ l)(L x),

where `f` is smooth, `g` and `h` are possibly nonsmooth and `l` is strongly
convex. Symbol `□` denotes the infimal convolution, and `L` is a linear mapping.
If `solver = VuCondat(args...)`, then the above problem is solved with

    solver(x0, y0; [f, g, h, l, L, beta_f, beta_l])

See documentation of `AFBA` for the list of keyword arguments.

References:

[1] Chambolle, Pock, "A First-Order Primal-Dual Algorithm for Convex Problems
with Applications to Imaging", Journal of Mathematical Imaging and Vision,
vol. 40, no. 1, pp. 120-145 (2011).

[2] Condat, "A primal–dual splitting method for convex optimization
involving Lipschitzian, proximable and linear composite terms",
Journal of Optimization Theory and Applications, vol. 158, no. 2,
pp 460-479 (2013).

[3] Vũ, "A splitting algorithm for dual monotone inclusions involving
cocoercive operators", Advances in Computational Mathematics, vol. 38, no. 3,
pp. 667-681 (2013).
"""
VuCondat(::Type{R}; kwargs...) where R = AFBA{R}(; kwargs..., theta=R(2.0))
VuCondat(; kwargs...) = VuCondat(Float64; kwargs...)