text:
mutable struct NSGA3 <: AbstractParameters
    N::Int
    η_cr::Float64
    p_cr::Float64
    η_m::Float64
    p_m::Float64
    partitions::Int
    reference_points::Array{Vector{Float64},1}
end

"""
    NSGA3(;
        N = 100,
        η_cr = 20,
        p_cr = 0.9,
        η_m = 20,
        p_m = 1.0 / D,
        partitions = 12,
        reference_points = Vector{Float64}[],
        information = Information(),
        options = Options(),
    )

Parameters for the metaheuristic NSGA-III.

Parameters:

- `N` Population size.
- `η_cr`  η for the crossover.
- `p_cr` Crossover probability.
- `η_m`  η for the mutation operator.
- `p_m` Mutation probability (1/D for D-dimensional problem by default).
- `reference_points` reference points usually generated by `gen_ref_dirs`.
- `partitions` number of Das and Dennis's reference points if `reference_points` is empty.

To use NSGA3, the output from the objective function should be a 3-touple
`(f::Vector, g::Vector, h::Vector)`, where `f` contains the objective functions,
`g` and `h` are inequality, equality constraints respectively.

A feasible solution is such that `g_i(x) ≤ 0 and h_j(x) = 0`.

```julia
using Metaheuristics

# Objective function, bounds, and the True Pareto front
f, bounds, pf = Metaheuristics.TestProblems.get_problem(:DTLZ2)

# define the parameters (use `NSGA3()` for using default parameters)
nsga3 = NSGA3(p_cr = 0.9)

# optimize
status = optimize(f, bounds, nsga3)

# show results
display(status)
```

"""
function NSGA3(;
    N = 100,
    η_cr = 20,
    p_cr = 0.9,
    η_m = 20,
    p_m = -1,
    partitions=12,
    reference_points=Vector{Float64}[],
    information = Information(),
    options = Options(),
)

    parameters = NSGA3(N, promote( Float64(η_cr), p_cr, η_m, p_m )..., partitions,
                      reference_points)
    Algorithm(
        parameters,
        information = information,
        options = options,
    )

end

function update_state!(
    status::State,
    parameters::NSGA3,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )

    I = randperm(parameters.N)
    for i = 1:parameters.N ÷ 2

        pa = status.population[I[2i-1]]
        pb = status.population[I[2i]]

        # crossover
        c1, c2 = SBX_crossover( get_position(pa), get_position(pb), problem.bounds,
                              parameters.η_cr, parameters.p_cr)

        # mutation
        polynomial_mutation!(c1,problem.bounds,parameters.η_m, parameters.p_m)
        polynomial_mutation!(c2,problem.bounds,parameters.η_m, parameters.p_m)

        # rapair solutions if necesary
        reset_to_violated_bounds!(c1, problem.bounds)
        reset_to_violated_bounds!(c2, problem.bounds)

        # evaluate children
        child1 = create_solution(c1, problem)

        child2 = create_solution(c2, problem)
        status.f_calls += 2

        # save children
        push!(status.population, child1)
        push!(status.population, child2)
    end

    # non-dominated sort, elitist removing via niching
    truncate_population_nsga3!(status.population,parameters.reference_points,parameters.N)

    stop_criteria!(status, parameters, problem, information, options)
end

function truncate_population_nsga3!(population, reference_points, N)
    fast_non_dominated_sort!(population)

    k = 1
    l = 1
    while l <= N
        k = population[l].rank
        l += 1
    end

    let k = k
        l = findlast(sol -> sol.rank == k, population)
        deleteat!(population, l+1:length(population))
    end

    l = findfirst(sol -> sol.rank == k, population)
    niching!(population, reference_points, N, l)
end

distance_point_to_rect(s, w) = @fastmath norm(s - (dot(w,s) / dot(w, w))*w  )

function associate!(nich, nich_freq, distance, F, reference_points, l)
    N = length(F)

    # find closest nich to corresponding solution
    for i = 1:N
        for j = 1:length(reference_points)
            d = distance_point_to_rect(F[i], reference_points[j])

            distance[i] < d && continue

            distance[i] = d
            nich[i] = j
        end

        # not associate last  front
        if i < l
            nich_freq[nich[i]] += 1
        end

    end
end

function hyperplane_normalization(population)
    M = length(fval(population[1]))

    ideal_point = ideal(population)
    nadir_point = fill(Inf, length(ideal_point))

    Fx = fvals(population) .- ideal_point'

    # identify extreme points
    extreme_points = zeros(Int, M)
    w = LinearAlgebra.I + fill(1e-6, M, M)

    for i in 1:M
        extreme_points[i] = argmin(nadir(Fx' ./ w[i,:]))
    end

    # check if intercepts can be obtained
    S = Fx[extreme_points,:]
    if LinearAlgebra.det(S) ≈ 0 # check if soluble matrix
        nadir_point = nadir(population)
    else
        hyperplane = S \ ones(M)
        intercepts = 1 ./ hyperplane # intercepts
        nadir_point = ideal_point + intercepts
    end

    ideal_point, nadir_point
end

function normalize(population)
    ideal_point, nadir_point = hyperplane_normalization(population)

    b = nadir_point - ideal_point

    # prevent division by zero
    mask = b .< eps()
    b[mask] .= eps()

    return [ (sol.f - ideal_point) ./ b for sol in population ]
end

# get_last_front(id, population) = findall(s -> s.rank == id, population)

pick_random(itr, item) = rand(findall(i -> i == item, itr))
find_item(itr, item)  = findall(i -> i == item, itr)

function niching!(population, reference_points, N, l)
    if length(population) == N
        return
    end

    F = normalize(population)
    k = l

    # allocate memory
    nich = zeros(Int, length(population))
    nich_freq = zeros(Int, length(reference_points))
    available_niches = ones(Bool, length(reference_points))
    distance = fill(Inf, length(population))

    # associate to niches
    associate!(nich, nich_freq, distance, F, reference_points, l)

    # keep last front
    last_front_id = k:length(population)#get_last_front(population[end].rank, population)
    last_front = population[last_front_id]
    deleteat!(population, last_front_id)

    # last front niches information
    niches_last_front = nich[last_front_id]
    distance_last_front = distance[last_front_id]

    # save survivors
    i = 1
    while k <= N
        mini = minimum(view(nich_freq, available_niches))
        # nich to be assigned
        j_hat = pick_random(nich_freq, mini)

        # candidate solutions
        I_j_hat = find_item( niches_last_front, j_hat )
        if isempty(I_j_hat)
            available_niches[j_hat] = false
            continue
        end

        if mini == 0
            ds = view(distance_last_front, I_j_hat)
            s = I_j_hat[argmin(ds)]
            push!(population, last_front[s])
        else
            s = rand(I_j_hat)
            push!(population, last_front[s])
        end

        nich_freq[j_hat] += 1
        deleteat!(last_front, s)
        deleteat!(niches_last_front, s)
        deleteat!(distance_last_front, s)

        k += 1

    end

    nothing
end

function initialize!(
    status,
    parameters::NSGA3,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
)
    D = size(problem.bounds, 2)

    if parameters.p_m < 0.0
        parameters.p_m = 1.0 / D
    end

    if options.iterations == 0
        options.iterations = 500
    end

    if options.f_calls_limit == 0
        options.f_calls_limit = options.iterations * parameters.N + 1
    end

    status = gen_initial_state(problem, parameters, information, options,status)

    if isempty(parameters.reference_points)
        options.debug && @info "Initializing reference points..."
        # number of objectives
        m = length(status.population[1].f)

        # initialize reference points
        parameters.reference_points = gen_weights(m, parameters.partitions)
    end

    status

end

function final_stage!(
    status::State,
    parameters::NSGA3,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )
    status.final_time = time()

end


text:
# complexspace.jl
#
# Defines the immutable ComplexSpace for a standard complex vector space C^n and is
# characterized by its dimension and whether or not it is the dual space. Tensors with
# ComplexSpace as index spaces make a distinction between covariant and contravariant indices.

# ComplexSpace:
#---------------
immutable ComplexSpace <: EuclideanSpace{ℂ}
  d::Int
  dual::Bool
end
ComplexSpace(d::Int) = ComplexSpace(d, false)
^(::Type{ℂ},d::Int) = ComplexSpace(d)

# Corresponding methods:
dim(V::ComplexSpace) = V.d
dual(V::ComplexSpace) = ComplexSpace(V.d, !V.dual)
cnumber(V::ComplexSpace) = ComplexSpace(1, V.dual)
cnumber(::Type{ComplexSpace}) = ComplexSpace(1)
iscnumber(V::ComplexSpace) = dim(V)==1

# Show methods
Base.show(io::IO, V::ComplexSpace) = print(io, V.dual ? "ℂ^$(V.d)*" : "ℂ^$(V.d)")

# direct sum of ComplexSpaces
directsum(V1::ComplexSpace, V2::ComplexSpace) = (V1.dual==V2.dual ? ComplexSpace(V1.d+V2.d, V1.dual) : throw(SpaceError("Direct sum of a vector space and its dual do not exist")))

# fusing and splitting ComplexSpaces
fuse(V1::ComplexSpace,V2::ComplexSpace,V::ComplexSpace) = dim(V1)*dim(V2)==dim(V)

# # basis and basisvector
# typealias ComplexBasisVector BasisVector{ComplexSpace,Int} # use integer from 1 to dim as identifier
# typealias ComplexBasis Basis{ComplexSpace}
#
# Base.length(B::ComplexBasis) = dim(space(B))
# Base.start(B::ComplexBasis) = 1
# Base.next(B::ComplexBasis, state::Int) = (EuclideanBasisVector(space(B),state),state+1)
# Base.done(B::ComplexBasis, state::Int) = state>length(B)
#
# Base.to_index(b::ComplexBasisVector) = b.identifier

ext: jl


text:
include("../optimizers/linesearch.jl")
include("../utils/functions.jl")
using Gadfly

function optimize_linear_constraints(f, A, f_g, f_h, yinit=0)
    """
    Convert constrained optimization problem min f(x) s.t. Ax = b
        to unconstrained problem min g(y) = f(B*y + x0), where x0 is a solution to Ax = b
        and B is a basis for the kernel of A.

    Defaults to using Newton's line search algorithm.

    f: constrained objective
    A: constraint matrix
    f_g: gradient of f
    f_h: hessian of f
    yinit: init y with a vector of all yinit values
    """

    m, n = size(A)
    b = zeros(m)
    x0 = A \ b
    B = nullspace(A)

    y0 = ones(size(B)[2]) .* yinit

    g(y) = f(B*y + x0)
    g_g(y) = B' * f_g(B*y + x0)
    g_h(y) = B' * f_h(B*y + x0) * B

    yvals = line_search(g, y0, g_g, g_h, "newton", 5000)
    xvals = [B*y + x0 for y in yvals]
    lambdas = [A' \ f_g(x) for x in xvals]

    return yvals, xvals, lambdas
end

function generate_constraint_matrix(n)
    A = ones(2,n)
    for i in 1:n
        if i % 2 == 1
            A[2,i] = -1
        end
    end
    return A
end

A = generate_constraint_matrix(30)
b = zeros(2)
x0 = A \ b; B = nullspace(A)
f = cute
f_g = cute_g
f_h = cute_h
g(y) = f(B*y + x0)
g_g(y) = B' * f_g(B*y + x0)
g_h(y) = B' * f_h(B*y + x0) * B

yvals, xvals, lambdas = optimize_linear_constraints(f, A, f_g, f_h, 0);

niters = length(xvals)
fx = [f(x) for x in xvals]
Gadfly.plot(x=1:niters, y=fx, Geom.line,
    Guide.xlabel("iteration"), Guide.ylabel("f(x)"), Guide.title("Value of function"),
    Scale.x_log10, Scale.y_log10)

grads = [norm(g_g(y),2) for y in yvals]
ratios = grads[2:niters,:]./grads[1:niters-1,:]
Gadfly.plot(x=1:niters-1, y=ratios, Geom.line,
Guide.xlabel("iteration"), Guide.ylabel("f_g(x)"), Guide.title("Gradient norm ratios"),
    Scale.x_log10, Scale.y_log10)

Gadfly.plot(x=1:niters, y=grads, Geom.line,
Guide.xlabel("iteration"), Guide.ylabel("f_g(x)"), Guide.title("Gradient norm"),
    Scale.x_log10, Scale.y_log10)

# Lagrange multipliers
lambdas[end]

A = generate_constraint_matrix(100)
x0 = A \ b; B = nullspace(A)
f = cute
f_g = cute_g
f_h = cute_h
g(y) = f(B*y + x0)
g_g(y) = B' * f_g(B*y + x0)
g_h(y) = B' * f_h(B*y + x0) * B

yvals, xvals, lambdas = optimize_linear_constraints(f, A, f_g, f_h, 0);
# Lagrange multipliers
println(lambdas[end])
# Final value
println(f(xvals[end]))

niters = length(xvals)
fx = [f(x) for x in xvals]
#gy = [g(y) for y in yvals]

Gadfly.plot(x=1:niters, y=fx, Geom.line,
    Guide.xlabel("iteration"), Guide.ylabel("f(x)"), Guide.title("Value of function"),
    Scale.x_log10, Scale.y_log10)

grads = [norm(g_g(y),2) for y in yvals]
ratios = grads[2:niters,:]./grads[1:niters-1,:]
Gadfly.plot(x=1:niters-1, y=ratios, Geom.line,
Guide.xlabel("iteration"), Guide.ylabel("f_g(x)"), Guide.title("Gradient norm ratios"),
    Scale.x_log10, Scale.y_log10)

Gadfly.plot(x=1:niters, y=grads, Geom.line,
Guide.xlabel("iteration"), Guide.ylabel("f_g(x)"), Guide.title("Gradient norm"),
    Scale.x_log10, Scale.y_log10)


text:
struct Span{P<:AbstractTracerProvider} <: AbstractSpan
    name::Ref{String}
    tracer::Tracer{P}
    span_context::SpanContext
    parent_span_context::Union{Nothing,SpanContext}
    kind::SpanKind
    start_time::UInt
    end_time::Ref{Union{Nothing,UInt}}
    attributes::DynamicAttrs
    links::Limited{Vector{Link}}
    events::Limited{Vector{Event}}
    status::Ref{SpanStatus}
end

function OpenTelemetryAPI.create_span(
    name::String,
    tracer::Tracer{<:TracerProvider};
    context = current_context(),
    kind = SPAN_KIND_INTERNAL,
    attributes = Dict{String,TAttrVal}(),
    links = Link[],
    events = OpenTelemetryAPI.Event[],
    start_time = UInt(time() * 10^9),
    is_remote = false,
)
    provider = tracer.provider
    parent_span_ctx = context |> current_span |> span_context
    if isnothing(parent_span_ctx)
        trace_id = generate_trace_id(provider.id_generator)
    else
        trace_id = parent_span_ctx.trace_id
    end

    attributes = DynamicAttrs(
        attributes;
        count_limit = provider.limit_info.span_attribute_count_limit,
        value_length_limit = provider.limit_info.span_attribute_value_length_limit,
    )

    links = Limited(links; limit = provider.limit_info.span_link_count_limit)
    events = Limited(events; limit = provider.limit_info.span_event_count_limit)

    sampling_result = should_sample(
        provider.sampler,
        context,
        trace_id,
        name,
        kind,
        attributes,
        links,
        isnothing(parent_span_ctx) ? TraceState() : parent_span_ctx.trace_state,
    )

    span_ctx = SpanContext(
        trace_id = trace_id,
        span_id = generate_span_id(provider.id_generator),
        is_remote = is_remote,
        trace_flag = TraceFlag(sampled = is_sampled(sampling_result)),
        trace_state = sampling_result.trace_state,
    )

    is_no_op_span = provider.is_shut_down[] || !is_recording(sampling_result)
    s = Span(
        Ref(name),
        tracer,
        is_no_op_span ? INVALID_SPAN_CONTEXT : span_ctx,
        parent_span_ctx,
        kind,
        start_time,
        Ref{Union{Nothing,UInt64}}(is_no_op_span ? start_time : nothing),
        attributes,
        links,
        events,
        Ref(SpanStatus(SPAN_STATUS_UNSET)),
    )
    if !is_no_op_span
        on_start!(tracer.provider.span_processor, s)
    end
    s
end

OpenTelemetryAPI.tracer(s::Span) = s.tracer
OpenTelemetryAPI.span_context(s::Span) = s.span_context
OpenTelemetryAPI.span_kind(s::Span) = s.kind
OpenTelemetryAPI.parent_span_context(s::Span) = s.parent_span_context
OpenTelemetryAPI.attributes(s::Span) = s.attributes
OpenTelemetryAPI.is_recording(s::Span) = isnothing(s.end_time[])
OpenTelemetryAPI.start_time(s::Span) = s.start_time
OpenTelemetryAPI.span_links(s::Span) = s.links
OpenTelemetryAPI.span_events(s::Span) = s.events
OpenTelemetryAPI.end_time(s::Span) = s.end_time[]

function OpenTelemetryAPI.end_span!(s::Span{<:TracerProvider}, t = UInt(time() * 10^9))
    if is_recording(s)
        s.end_time[] = t
        on_end!(s.tracer.provider.span_processor, s)
    else
        @warn "the span is not recording."
    end
end

function OpenTelemetryAPI.span_status!(s::Span, code::SpanStatusCode, description = nothing)
    if is_recording(s)
        if s.status[].code === SPAN_STATUS_OK
            # no further updates
        else
            if code === SPAN_STATUS_UNSET
                # ignore
            else
                s.status[] = SpanStatus(code, description)
            end
        end
    else
        @warn "the span is not recording."
    end
end

OpenTelemetryAPI.span_status(s::Span) = s.status[]

OpenTelemetryAPI.span_name(s::Span) = s.name[]

function OpenTelemetryAPI.span_name!(s::Span, name::String)
    if is_recording(s)
        s.name[] = name
    else
        @warn "the span is not recording."
    end
end

function Base.push!(s::Span, ex::Exception; is_rethrow_followed = false)
    msg_io = IOBuffer()
    showerror(msg_io, ex)
    msg = String(take!(msg_io))

    st_io = IOBuffer()
    showerror(st_io, CapturedException(ex, catch_backtrace()))
    st = String(take!(st_io))

    attrs = StaticAttrs((;
        Symbol("exception.type") => string(typeof(ex)),
        Symbol("exception.type") => string(typeof(ex)),
        Symbol("exception.message") => msg,
        Symbol("exception.stacktrace") => st,
        Symbol("exception.escaped") => is_rethrow_followed,
    ))

    push!(s, Event(name = "exception", attributes = attrs))
end


text:
export PlainProbCircuit,
    PlainProbLeafNode, PlainProbInnerNode,
    PlainProbLiteralNode, PlainMulNode, PlainSumNode

#####################
# Plain probabilistic circuit nodes
#####################

"Root of the plain probabilistic circuit node hierarchy"
abstract type PlainProbCircuit <: ProbCircuit end

"A probabilistic leaf node"
abstract type PlainProbLeafNode <: PlainProbCircuit end

"A probabilistic inner node"
abstract type PlainProbInnerNode <: PlainProbCircuit end

"A probabilistic literal node"
mutable struct PlainProbLiteralNode <: PlainProbLeafNode
    literal::Lit
    data
    counter::UInt32
    PlainProbLiteralNode(l) = new(l, nothing, 0)
end

"A probabilistic conjunction node (multiplication node)"
mutable struct PlainMulNode <: PlainProbInnerNode
    children::Vector{PlainProbCircuit}
    data
    counter::UInt32
    PlainMulNode(children) = begin
        new(convert(Vector{PlainProbCircuit}, children), nothing, 0)
    end
end

"A probabilistic disjunction node (summation node)"
mutable struct PlainSumNode <: PlainProbInnerNode
    children::Vector{PlainProbCircuit}
    log_probs::Vector{Float64}
    data
    counter::UInt32
    PlainSumNode(c) = begin
        new(c, log.(ones(Float64, length(c)) / length(c)), nothing, 0)
    end
end

#####################
# traits
#####################

import LogicCircuits.GateType # make available for extension

@inline GateType(::Type{<:PlainProbLiteralNode}) = LiteralGate()
@inline GateType(::Type{<:PlainMulNode}) = ⋀Gate()
@inline GateType(::Type{<:PlainSumNode}) = ⋁Gate()

#####################
# methods
#####################

import LogicCircuits: children # make available for extension
@inline children(n::PlainProbInnerNode) = n.children

"Count the number of parameters in the node"
@inline num_parameters_node(n::PlainSumNode) = num_children(n)

#####################
# constructors and conversions
#####################

function multiply(arguments::Vector{<:PlainProbCircuit};
                 reuse=nothing)
    @assert length(arguments) > 0
    reuse isa PlainMulNode && children(reuse) == arguments && return reuse
    return PlainMulNode(arguments)
end

function summate(arguments::Vector{<:PlainProbCircuit};
                    reuse=nothing)
    @assert length(arguments) > 0
    reuse isa PlainSumNode && children(reuse) == arguments && return reuse
    return PlainSumNode(arguments)
end

# claim `PlainProbCircuit` as the default `ProbCircuit` implementation
compile(::Type{ProbCircuit}, args...) =
    compile(PlainProbCircuit, args...)

compile(::Type{<:PlainProbCircuit}, l::Lit) =
    PlainProbLiteralNode(l)

function compile(::Type{<:PlainProbCircuit}, circuit::LogicCircuit)
    f_con(n) = error("Cannot construct a probabilistic circuit from constant leafs: first smooth and remove unsatisfiable branches.")
    f_lit(n) = compile(PlainProbCircuit, literal(n))
    f_a(_, cns) = multiply(cns)
    f_o(_, cns) = summate(cns)
    foldup_aggregate(circuit, f_con, f_lit, f_a, f_o, PlainProbCircuit)
end

import LogicCircuits: fully_factorized_circuit #extend

fully_factorized_circuit(::Type{ProbCircuit}, n::Int) =
    fully_factorized_circuit(PlainProbCircuit, n)

function fully_factorized_circuit(::Type{<:PlainProbCircuit}, n::Int)
    ff_logic_circuit = fully_factorized_circuit(PlainLogicCircuit, n)
    compile(PlainProbCircuit, ff_logic_circuit)
end


text:
function get_incompatible_devices(devices_template::Dict{Symbol, DeviceModel})
    incompatible_device_types = Vector{DataType}()
    for model in values(devices_template)
        formulation = get_formulation(model)
        if formulation == FixedOutput
            if !isempty(get_services(model))
                @info "$(formulation) for $(get_device_type(model)) is not compatible with the provision of reserve services"
            end
            push!(incompatible_device_types, get_device_type(model))
        end
    end
    return incompatible_device_types
end

function construct_services!(
    psi_container::PSIContainer,
    sys::PSY.System,
    services_template::Dict{Symbol, ServiceModel},
    devices_template::Dict{Symbol, DeviceModel},
)
    isempty(services_template) && return
    incompatible_device_types = get_incompatible_devices(devices_template)

    function _construct_valid_services!(service_model::ServiceModel)
        @debug "Building $(service_model.service_type) with $(service_model.formulation) formulation"
        services = service_model.service_type[]
        if validate_services!(
            service_model.service_type,
            services,
            incompatible_device_types,
            sys,
        )
            construct_service!(
                psi_container,
                services,
                sys,
                service_model,
                devices_template,
                incompatible_device_types,
            )
        end
    end

    groupservice = nothing

    for (key, service_model) in services_template
        if service_model.formulation === GroupReserve  # group service needs to be constructed last
            groupservice = key
            continue
        end
        _construct_valid_services!(service_model)
    end
    groupservice === nothing || _construct_valid_services!(services_template[groupservice])
    return
end

function construct_service!(
    psi_container::PSIContainer,
    services::Vector{SR},
    sys::PSY.System,
    model::ServiceModel{SR, RangeReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Vector{<:DataType},
) where {SR <: PSY.Reserve}
    services_mapping = PSY.get_contributing_device_mapping(sys)
    time_steps = model_time_steps(psi_container)
    names = [PSY.get_name(s) for s in services]

    if model_has_parameters(psi_container)
        container = add_param_container!(
            psi_container,
            UpdateRef{SR}("service_requirement", "get_requirement"),
            names,
            time_steps,
        )
    end

    add_cons_container!(
        psi_container,
        make_constraint_name(REQUIREMENT, SR),
        names,
        time_steps,
    )

    for service in services
        contributing_devices =
            services_mapping[(type = SR, name = PSY.get_name(service))].contributing_devices
        if !isempty(incompatible_device_types)
            contributing_devices =
                [d for d in contributing_devices if typeof(d) ∉ incompatible_device_types]
        end
        #Services without contributing devices should have been filtered out in the validation
        @assert !isempty(contributing_devices)
        #Variables
        add_variables!(psi_container, ActiveServiceVariable, service, contributing_devices)
        # Constraints
        service_requirement_constraint!(psi_container, service, model)
        modify_device_model!(devices_template, model, contributing_devices)

        # Cost Function
        cost_function!(psi_container, service, model)
    end
    return
end

function construct_service!(
    psi_container::PSIContainer,
    services::Vector{SR},
    sys::PSY.System,
    model::ServiceModel{SR, StepwiseCostReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Vector{<:DataType},
) where {SR <: PSY.Reserve}
    services_mapping = PSY.get_contributing_device_mapping(sys)
    time_steps = model_time_steps(psi_container)
    names = [PSY.get_name(s) for s in services]
    add_variables!(psi_container, ServiceRequirementVariable, services)
    add_cons_container!(
        psi_container,
        make_constraint_name(REQUIREMENT, SR),
        names,
        time_steps,
    )

    for service in services
        contributing_devices =
            services_mapping[(
                type = typeof(service),
                name = PSY.get_name(service),
            )].contributing_devices
        if !isempty(incompatible_device_types)
            contributing_devices =
                [d for d in contributing_devices if typeof(d) ∉ incompatible_device_types]
        end
        #Variables
        add_variables!(psi_container, ActiveServiceVariable, service, contributing_devices)
        # Constraints
        service_requirement_constraint!(psi_container, service, model)
        modify_device_model!(devices_template, model, contributing_devices)

        # Cost Function
        cost_function!(psi_container, service, model.formulation)
    end
    return
end

function construct_service!(
    psi_container::PSIContainer,
    services::Vector{PSY.AGC},
    sys::PSY.System,
    ::ServiceModel{PSY.AGC, T},
    devices_template::Dict{Symbol, DeviceModel},
    ::Vector{<:DataType},
) where {T <: AbstractAGCFormulation}
    #Order is important in the addition of these variables
    for device_model in devices_template
        #TODO: make a check for the devices' models
    end
    agc_areas = [PSY.get_area(agc) for agc in services]
    areas = PSY.get_components(PSY.Area, sys)
    for area in areas
        if area ∉ agc_areas
            #    throw(IS.ConflictingInputsError("All area most have an AGC service assigned in order to model the System's Frequency regulation"))
        end
    end
    add_variables!(psi_container, SteadyStateFrequencyDeviation)
    add_variables!(psi_container, AreaMismatchVariable, areas)
    add_variables!(psi_container, SmoothACE, areas)
    add_variables!(psi_container, LiftVariable, areas)
    add_variables!(psi_container, ActivePowerVariable, areas)
    add_variables!(psi_container, DeltaActivePowerUpVariable, areas)
    add_variables!(psi_container, DeltaActivePowerDownVariable, areas)
    #add_variables!(psi_container, AdditionalDeltaActivePowerUpVariable, areas)
    #add_variables!(psi_container, AdditionalDeltaActivePowerDownVariable, areas)
    balancing_auxiliary_variables!(psi_container, sys)

    absolute_value_lift(psi_container, areas)
    frequency_response_constraint!(psi_container, sys)
    area_control_init(psi_container, services)
    smooth_ace_pid!(psi_container, services)
    aux_constraints!(psi_container, sys)
end

"""
    Constructs a service for StaticReserveGroup.
"""
function construct_service!(
    psi_container::PSIContainer,
    services::Vector{SR},
    ::PSY.System,
    model::ServiceModel{SR, GroupReserve},
    ::Dict{Symbol, DeviceModel},
    ::Vector{<:DataType},
) where {SR <: PSY.StaticReserveGroup}
    time_steps = model_time_steps(psi_container)
    names = [PSY.get_name(s) for s in services]

    if model_has_parameters(psi_container)
        container = add_param_container!(
            psi_container,
            UpdateRef{SR}("service_requirement", "get_requirement"),
            names,
            time_steps,
        )
    end

    add_cons_container!(
        psi_container,
        make_constraint_name(REQUIREMENT, SR),
        names,
        time_steps,
    )

    for service in services
        contributing_services = PSY.get_contributing_services(service)

        # check if variables exist
        check_activeservice_variables(psi_container, contributing_services)
        # Constraints
        service_requirement_constraint!(
            psi_container,
            service,
            model,
            contributing_services,
        )
    end
    return
end


text:
module GVFSRHordes
import Lazy: @forward

export SFHorde, SRHorde
using GVFHordes
using ..Curiosity: FeatureCreator

struct SFHorde{T<:AbstractGVF} <: GVFHordes.AbstractHorde
    gvfs::Vector{T}
end

# mutable struct SRHorde{T<:GVFHordes.AbstractHorde} <: GVFHordes.AbstractHorde
mutable struct SRHorde <: GVFHordes.AbstractHorde
    PredHorde::GVFHordes.AbstractHorde
    SFHorde::SFHorde
    num_tasks::Int
    num_SFs::Int
    # state_constructor::Union{AbstractFeatureProjector, Curiosity.TileCoder}
    state_constructor::FeatureCreator
    # state_constructor::Function
    function SRHorde(prediction_horde::GVFHordes.AbstractHorde, successor_feature_horde::SFHorde, num_SFs, state_constructor)
        new(prediction_horde, successor_feature_horde, length(prediction_horde), num_SFs, state_constructor)
    end
end

function assign_prediction_horde!(h::SRHorde, new_pred_horde)
    SRHorde.PredHorde = new_pred_horde
end

function Base.length(h::SRHorde)
    return length(h.PredHorde) + length(h.SFHorde)
end

# function Base.get(gvfh::SRHorde, state_t, action_t, state_tp1, action_tp1, preds_tp1)
function Base.get(gvfh::SRHorde; kwargs...)
    C, discounts, pi = get(gvfh.PredHorde; kwargs...)
    # The reward discounts should always be 0 as they used for supervised learning prediction.
    @assert sum(discounts) == 0
    # For efficiency, the discount and pi for a given SF task should all have the same values. Therefore, to
    # improve speed, a single call to an GVF within each block should be sufficient.
    state_action_feature_length = Int( length(gvfh.SFHorde) / gvfh.num_SFs)
    inds_per_task = collect(1:state_action_feature_length:length(gvfh.SFHorde))

    constructed_state = gvfh.state_constructor(kwargs[:state_t], kwargs[:action_t], kwargs[:state_tp1])

    # @show constructed_state, kwargs[:action_t]
    #NOTE: Commenting out potential optimization as it is hacky and not too big of saving.
    # nz_inds = findall(!iszero, constructed_state)
    # num_actions = 4
    # C_SF_nz_ind = (nz_inds .- 1) * num_actions .+ kwargs[:action_t]
    # C_SF_per_SF = zeros(state_action_feature_length)
    # C_SF_per_SF[C_SF_nz_ind] .= 1
    # C_SF_fast = repeat(C_SF_per_SF, gvfh.num_SFs)
    #
    # cumulant_inds_per_task = state_action_feature_length
    # # C_SF = C_SF_fast
    C_SF = map(gvf -> get(cumulant(gvf); constructed_state_t = constructed_state, kwargs...), gvfh.SFHorde.gvfs)
    # #
    # if C_SF != C_SF_fast
    #     throw(ArgumentError("Fast version is not the same"))
    #     @show size(C_SF), size(C_SF_fast)
    # end
    # C_SF = zeros(length(gvfh.SFHorde))
    #
    # @show findall(!iszero, C_SF)
    # println()
    # println()

    # unique_discounts_SF = map(gvf -> get(discount(gvf), state_t, action_t, state_tp1, action_tp1, preds_tp1), gvfh.SFHorde.gvfs[inds_per_task])
    unique_discounts_SF = map(gvf -> get(discount(gvf); kwargs...), gvfh.SFHorde.gvfs[inds_per_task])
    # unique_pi_SF = map(gvf -> get(policy(gvf), state_t, action_t), gvfh.SFHorde.gvfs[inds_per_task])
    unique_pi_SF = map(gvf -> get(policy(gvf); kwargs...), gvfh.SFHorde.gvfs[inds_per_task])

    discounts_SF = repeat(unique_discounts_SF, inner = state_action_feature_length)
    pi_SF = repeat(unique_pi_SF, inner = state_action_feature_length)

    return vcat(C,C_SF), vcat(discounts,discounts_SF), vcat(pi,pi_SF)
end

# Base.get(gvfh::SRHorde, state_tp1, preds_tp1) =
#     get(gvfh::SRHorde, nothing, nothing, state_tp1, nothing, preds_tp1)
#
# Base.get(gvfh::SRHorde, state_t, action_t, state_tp1) =
#     get(gvfh::SRHorde, state_t, action_t, state_tp1, nothing, nothing)
#
# Base.get(gvfh::SRHorde, state_t, action_t, state_tp1, preds_tp1) =
#     get(gvfh::SRHorde, state_t, action_t, state_tp1, nothing, preds_tp1)

Base.get(gvfh::SRHorde, state_tp1, preds_tp1) =
    get(gvfh::SRHorde; state_t = nothing,
    action_t = nothing,
    state_tp1 = state_tp1,
    action_tp1 = nothing,
    preds_tp1 = nothing)

Base.get(gvfh::SRHorde, state_t, action_t, state_tp1) =
    get(gvfh::SRHorde;
    state_t = state_t,
    action_t = action_t,
    state_tp1 = state_tp1,
    action_tp1 = nothing,
    preds_tp1 = nothing)

Base.get(gvfh::SRHorde, state_t, action_t, state_tp1, preds_tp1) =
    get(gvfh::SRHorde;
    state_t = state_t,
    action_t = action_t,
    state_tp1 = state_tp1,
    action_tp1 = nothing,
    preds_tp1 = preds_tp1)

merge(h1::SFHorde, h2::SFHorde) =
    SFHorde([deepcopy(h1.gvfs); deepcopy(h2.gvfs)])

function cumulant(gvf::AbstractGVF) end
function discount(gvf::AbstractGVF) end
function policy(gvf::AbstractGVF) end

cumulant(gvf::GVF) = gvf.cumulant
discount(gvf::GVF) = gvf.discount
policy(gvf::GVF) = gvf.policy

function Base.get(gvfh::SFHorde, state_t, action_t, state_tp1, action_tp1, preds_tp1)
    C = map(gvf -> get(cumulant(gvf), state_t, action_t, preds_tp1), gvfh.gvfs)
    Γ = map(gvf -> get(discount(gvf), state_t, action_t, state_tp1, action_tp1, preds_tp1), gvfh.gvfs)
    Π_probs = map(gvf -> get(policy(gvf), state_t, action_t), gvfh.gvfs)
    return C, Γ, Π_probs
end

function Base.get!(C::Array{T, 1}, Γ::Array{F, 1}, Π_probs::Array{H, 1}, gvfh::SFHorde, state_t, action_t, state_tp1, action_tp1, preds_tp1) where {T, F, H}
    C .= map(gvf -> get(cumulant(gvf), state_t, action_t, preds_t), gvfh.gvfs)
    Γ .= map(gvf -> get(discount(gvf), state_t, action_t, state_tp1, action_tp1, preds_tp1), gvfh.gvfs)
    Π_probs .= map(gvf -> get(policy(gvf), state_t, action_t), gvfh.gvfs)
    return C, Γ, Π_probs
end

Base.get(gvfh::SFHorde, state_tp1, preds_tp1) =
    get(gvfh::SFHorde; state_t = nothing,
    action_t = nothing,
    state_tp1 = state_tp1,
    action_tp1 = nothing,
    preds_tp1 = nothing)

Base.get(gvfh::SFHorde, state_t, action_t, state_tp1) =
    get(gvfh::SFHorde;
    state_t = state_t,
    action_t = action_t,
    state_tp1 = state_tp1,
    action_tp1 = nothing,
    preds_tp1 = nothing)

Base.get(gvfh::SFHorde, state_t, action_t, state_tp1, preds_tp1) =
    get(gvfh::SFHorde;
    state_t = state_t,
    action_t = action_t,
    state_tp1 = state_tp1,
    action_tp1 = nothing,
    preds_tp1 = preds_tp1)
#
# Base.get(gvfh::SFHorde, state_tp1, preds_tp1) =
#     get(gvfh::SFHorde, nothing, nothing, state_tp1, nothing, preds_tp1)
#
# Base.get(gvfh::SFHorde, state_t, action_t, state_tp1) =
#     get(gvfh::SFHorde, state_t, action_t, state_tp1, nothing, nothing)
#
# Base.get(gvfh::SFHorde, state_t, action_t, state_tp1, preds_tp1) =
#     get(gvfh::SFHorde, state_t, action_t, state_tp1, nothing, preds_tp1)

@forward SFHorde.gvfs Base.length
end #end of module

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using Statistics
using JLD2
using Printf
using GLMakie
using CUDA
using Oceananigans
using Oceananigans.Coriolis: HydrostaticSphericalCoriolis
using Oceananigans.Architectures: arch_array
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
using Oceananigans.TurbulenceClosures: HorizontallyCurvilinearAnisotropicDiffusivity
using Oceananigans.Units
using Oceananigans.Operators: Δzᵃᵃᶜ

#####
##### Grid
#####

latitude = (-84.375, 84.375)
Δφ = latitude[2] - latitude[1]

# 2.8125 degree resolution
Nx = 128
Ny = 60
Nz = 18

arch = CPU()
reference_density = 1035

#####
##### Load forcing files roughly from CORE2 paper
##### (Probably https://data1.gfdl.noaa.gov/nomads/forms/core/COREv2.html)
#####

bathymetry_path = "bathy_128x60var4.bin"
east_west_stress_path = "off_TAUXvar1.bin"
north_south_stress_path = "off_TAUY.bin"
sea_surface_temperature_path="sst25_128x60x12.bin"

bytes = sizeof(Float32) * Nx * Ny
bathymetry = reshape(bswap.(reinterpret(Float32, read(bathymetry_path, bytes))), (Nx, Ny))
τˣ = - reshape(bswap.(reinterpret(Float32, read(east_west_stress_path, bytes))), (Nx, Ny)) ./ reference_density
τʸ = - reshape(bswap.(reinterpret(Float32, read(north_south_stress_path, bytes))), (Nx, Ny)) ./ reference_density
target_sea_surface_temperature = reshape(bswap.(reinterpret(Float32, read(sea_surface_temperature_path, bytes))), (Nx, Ny))

bathymetry = arch_array(arch, bathymetry)
τˣ = arch_array(arch, τˣ)
τʸ = arch_array(arch, τʸ)
target_sea_surface_temperature = T★ = arch_array(arch, target_sea_surface_temperature)

H = 3600.0
# bathymetry = - H .* (bathymetry .< -10)
# H = - minimum(bathymetry)

# A spherical domain
@show underlying_grid = LatitudeLongitudeGrid(arch,
                                              size = (Nx, Ny, Nz),
                                              longitude = (-180, 180),
                                              latitude = latitude,
                                              halo = (3, 3, 3),
                                              z = (-H, 0))

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry))

#####
##### Physics and model setup
#####

νh = 1e+5
νz = 1e-2
κh = 1e+3
κz = 1e-4

background_diffusivity = HorizontallyCurvilinearAnisotropicDiffusivity(νh=νh, νz=νz, κh=κh, κz=κz)
convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0)

#####
##### Boundary conditions / constant-in-time surface forcing
#####

Δz_top = CUDA.@allowscalar Δzᵃᵃᶜ(1, 1, grid.Nz, grid.grid)
Δz_bottom = CUDA.@allowscalar Δzᵃᵃᶜ(1, 1, 1, grid.grid)

@inline surface_temperature_relaxation(i, j, grid, clock, fields, p) = @inbounds p.λ * (fields.T[i, j, grid.Nz] - p.T★[i, j])

T_surface_relaxation_bc = FluxBoundaryCondition(surface_temperature_relaxation,
                                                discrete_form = true,
                                                parameters = (λ = Δz_top/3days, T★ = target_sea_surface_temperature))

u_wind_stress_bc = FluxBoundaryCondition(τˣ)
v_wind_stress_bc = FluxBoundaryCondition(τʸ)

@inline u_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, 1]
@inline v_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, 1]

# Linear bottom drag:
μ = 1 / 10days * Δz_bottom

u_bottom_drag_bc = FluxBoundaryCondition(u_bottom_drag, discrete_form = true, parameters = μ)
v_bottom_drag_bc = FluxBoundaryCondition(v_bottom_drag, discrete_form = true, parameters = μ)

u_bcs = FieldBoundaryConditions(top = u_wind_stress_bc, bottom = u_bottom_drag_bc)
v_bcs = FieldBoundaryConditions(top = v_wind_stress_bc, bottom = v_bottom_drag_bc)
T_bcs = FieldBoundaryConditions(top = T_surface_relaxation_bc)

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    architecture = arch,
                                    #free_surface = ImplicitFreeSurface(maximum_iterations=10),
                                    #free_surface = ImplicitFreeSurface(),
                                    momentum_advection = VectorInvariant(),
                                    tracer_advection = WENO5(),
                                    coriolis = HydrostaticSphericalCoriolis(),
                                    boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs),
                                    buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=2e-4, β=0.0)),
                                    tracers = (:T, :S),
                                    closure = (background_diffusivity, convective_adjustment))

#=
#####
##### Visualize the setup
#####

setup_fig = Figure(resolution = (1200, 900))

free_surface = model.free_surface
∫Ax = free_surface.implicit_step_solver.vertically_integrated_lateral_areas.xᶠᶜᶜ
∫Ay = free_surface.implicit_step_solver.vertically_integrated_lateral_areas.yᶜᶠᶜ

h_cpu = Array(bathymetry)
τˣ_cpu = Array(τˣ)
τʸ_cpu = Array(τʸ)
T★_cpu = Array(T★)
∫Ax_cpu = Array(parent(∫Ax))[:, :, 1]
∫Ay_cpu = Array(parent(∫Ay))[:, :, 1]

ax_b = Axis(setup_fig[1, 1], title="Bathymetry (m)")
hm_b = heatmap!(ax_b, h_cpu)
cb_b = Colorbar(setup_fig[1, 2], hm_b)

ax_T★ = Axis(setup_fig[1, 3], title="Target sea surface temperature (ᵒC)")
hm_T★ = heatmap!(ax_T★, T★_cpu)
cb_T★ = Colorbar(setup_fig[1, 4], hm_T★)

ax_x = Axis(setup_fig[2, 1], title="Vertically-integrated x-face area (m³)")
hm_x = heatmap!(ax_x, ∫Ax_cpu)
cb_x = Colorbar(setup_fig[2, 2], hm_x)

ax_y = Axis(setup_fig[2, 3], title="Vertically-integrated y-face area (m³)")
hm_y = heatmap!(ax_y, ∫Ay_cpu)
cb_y = Colorbar(setup_fig[2, 4], hm_y)

ax_τˣ = Axis(setup_fig[3, 1], title="East-west wind stress (m² s⁻²)")
hm_τˣ = heatmap!(ax_τˣ, τˣ_cpu)
cb_τˣ = Colorbar(setup_fig[3, 2], hm_τˣ)

ax_τʸ = Axis(setup_fig[3, 3], title="North-south wind stress (m² s⁻²)")
hm_τʸ = heatmap!(ax_τʸ, τʸ_cpu)
cb_τʸ = Colorbar(setup_fig[3, 4], hm_τʸ)

display(setup_fig)
=#

#####
##### Initial condition:
#####

u, v, w = model.velocities
η = model.free_surface.η

T = model.tracers.T
T .= 5

S = model.tracers.S
S .= 30

#####
##### Simulation setup
#####

# Time-scale for gravity wave propagation across the smallest grid cell
g = model.free_surface.gravitational_acceleration
gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed

#minimum_Δx = abs(grid.radius * cosd(maximum(abs, grid.φᵃᶜᵃ[1:grid.Ny])) * deg2rad(grid.Δλ))
#minimum_Δy = abs(grid.radius * deg2rad(grid.Δφ))
#wave_propagation_time_scale = min(minimum_Δx, minimum_Δy) / gravity_wave_speed

if model.free_surface isa ExplicitFreeSurface
    Δt = 60seconds #0.2 * minimum_Δx / gravity_wave_speed
else
    Δt = 20minutes
end

simulation = Simulation(model, Δt = Δt, stop_time = 1day)

start_time = [time_ns()]

function progress(sim)
    wall_time = (time_ns() - start_time[1]) * 1e-9

    η = model.free_surface.η

    if model.free_surface isa ExplicitFreeSurface
        @info @sprintf("Time: % 12s, iteration: %d, max(|η|): %.2e m, wall time: %s",
                       prettytime(sim.model.clock.time),
                       sim.model.clock.iteration,
                       maximum(abs, η),
                       prettytime(wall_time))
    else
        @info @sprintf("Time: % 12s, iteration: %d, free surface iterations: %d, max(|η|): %.2e m, wall time: %s",
                       prettytime(sim.model.clock.time),
                       sim.model.clock.iteration,
                       sim.model.free_surface.implicit_step_solver.preconditioned_conjugate_gradient_solver.iteration,
                       maximum(abs, η),
                       prettytime(wall_time))
    end

    start_time[1] = time_ns()

    return nothing
end

simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

output_fields = merge(model.velocities, model.tracers, (; η=model.free_surface.η))
output_prefix = "global_lat_lon_$(grid.Nx)_$(grid.Ny)_$(grid.Nz)"

simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields,
                                                      schedule = TimeInterval(10days),
                                                      prefix = output_prefix,
                                                      force = true)

# Let's goo!
@info "Running with Δt = $(prettytime(simulation.Δt))"

run!(simulation)

@info """

    Background diffusivity: $background_diffusivity
    Free surface: $(typeof(model.free_surface).name.wrapper)
    Time step: $(prettytime(Δt))

"""
#   Simulation took $(prettytime(simulation.run_time))
#   Minimum wave propagation time scale: $(prettytime(wave_propagation_time_scale))

#####
##### Visualize solution
#####

η_cpu = Array(interior(model.free_surface.η))[:, :, 1]
T_cpu = Array(interior(model.tracers.T))[:, :, end]
u_cpu = Array(interior(model.velocities.u))[:, :, end]
v_cpu = Array(interior(model.velocities.v))[:, :, end]

max_η = maximum(abs, η_cpu)
max_u = maximum(abs, u_cpu)
max_v = maximum(abs, v_cpu)

max_T = maximum(T_cpu)
min_T = minimum(T_cpu)

solution_fig = Figure(resolution = (1200, 600))

ax_η = Axis(solution_fig[1, 1], title="Free surface displacement (m)")
hm_η = heatmap!(ax_η, η_cpu, colorrange=(-max_η, max_η), colormap=:balance)
cb_η = Colorbar(solution_fig[1, 2], hm_η)

ax_T = Axis(solution_fig[2, 1], title="Sea surface temperature (ᵒC)")
hm_T = heatmap!(ax_T, T_cpu, colorrange=(min_T, max_T), colormap=:thermal)
cb_T = Colorbar(solution_fig[2, 2], hm_T)

ax_u = Axis(solution_fig[1, 3], title="East-West velocity (m s⁻¹)")
hm_u = heatmap!(ax_u, u_cpu, colorrange=(-max_u, max_u), colormap=:balance)
cb_u = Colorbar(solution_fig[1, 4], hm_u)

ax_v = Axis(solution_fig[2, 3], title="North-South velocity (m s⁻¹)")
hm_v = heatmap!(ax_v, v_cpu, colorrange=(-max_v, max_v), colormap=:balance)
cb_v = Colorbar(solution_fig[2, 4], hm_v)

display(solution_fig)


text:
module Sudoku

# package code goes here

using JuMP, Cbc, ChooseOptimizer

export sudoku, sudoku_check

"""
`sudoku(A::Matrix)` solves the Sudoku puzzle given by the matrix `A`.
Here `A` is a 9-by-9 integer matrix whose nonzero entries are the
given values of a Sudoku puzzle and whose zero values are the blanks.
"""
function sudoku(A::Matrix{Int})::Matrix{Int}
    nn = size(A)[1] # number of items per row/column/block
    n = Int(sqrt(nn)) # size of the problem: normal sudoku n = 3
    # MOD = Model(JuMP.optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))
    set_solver_verbose(false)
    MOD = Model(get_solver())

    # variable X[i,j,k] = 1 means there's a k in cell (i,j)
    @variable(MOD, X[1:nn, 1:nn, 1:nn], Bin)

    ########### model constraints ##################
    @constraints(
        MOD,
        begin
            oneKPerCell, sum(X[1:nn, 1:nn, k] for k ∈ 1:nn) .== 1
            oneKPerRow, sum(X[1:nn, j, 1:nn] for j ∈ 1:nn) .== 1
            oneKPerCol, sum(X[i, 1:nn, 1:nn] for i ∈ 1:nn) .== 1
            oneKPerBlock,
            sum.(X[i:i+n-1, j:j+n-1, k] for k ∈ 1:nn, i ∈ 1:n:nn, j ∈ 1:n:nn) .== 1
        end
    )

    # Process values in A
    nonZeroIndices = findall(i -> i != 0, A)
    for i in nonZeroIndices
        @constraint(MOD, X[i, A[i]] == 1)
    end

    # now solve and extract the solution
    optimize!(MOD)
    status = Int(termination_status(MOD))
    if status != 1
        error("No solution to this Sudoku puzzle")
    end

    return Int.(sum(value.(X[:, :, k]) * k for k = 1:nn))
end

"""
`sudoku_check(A)` checks if a 9x9 matrix is a valid solution to a Sudoku
puzzle. That is, we check that every row, column, and 3x3 submatrix contains
the values 1 through 9 exactly once each.
"""
function sudoku_check(A::Matrix{Int})::Bool
    nn = size(A)[1] # number of items per row/column/block
    n = Int(sqrt(nn)) # size of the problem: normal sudoku n = 3

    # All items in a row are 1:9
    if !all(mapslices(aux -> sort(aux) == 1:nn, A, dims = 1))
        return false
    end

    # All items in a column are 1:9
    if !all(mapslices(aux -> sort(aux) == 1:nn, A, dims = 2))
        return false
    end

    # Finally we check each block, if we get this far, the sucoku will be correct
    # or not decided by this condition
    return all(
        aux -> sort(aux) == 1:nn,
        vec.([A[i:i+n-1, j:j+n-1] for i ∈ 1:n:nn, j ∈ 1:n:nn]),
    )
end

include("sudoku_print.jl")
include("examples.jl")

end # module


text:
module Tokens

abstract type AbstractToken end

struct Num <: AbstractToken
    value::Number
end

struct RBrace <: AbstractToken
end

struct LBrace <: AbstractToken
end

abstract type AbstractOperator <: AbstractToken end

struct Plus <: AbstractOperator
end

struct Minus <: AbstractOperator
end

struct Mul <: AbstractOperator
end

struct Div <: AbstractOperator
end

struct Pow <: AbstractOperator
end

end

max_stars_repo_name: MhmCats/AsciiMath

max_issues_repo_name: MhmCats/AsciiMath

max_forks_repo_name: MhmCats/AsciiMath


text:
"""
    expand_poly(x, [degree = 5], [obsdim]) -> Matrix

Perform a simple polynomial basis expansion of the given `degree`
for the vector `x`.

The optional parameter `obsdim` can be used to specify which
dimension of the resulting `Matrix` should correspond to the
observations. By default each column will denote an observation,
which means that the resulting `Matrix` would have the size
`(degree, length(x))`.

```@jldoctest
julia> expand_poly(1:5, degree = 3)
3×5 Array{Int64,2}:
 1  2   3   4    5
 1  4   9  16   25
 1  8  27  64  125

julia> expand_poly(1f0:5f0, 3) # positional arguments are type-stable
3×5 Array{Float32,2}:
 1.0  2.0   3.0   4.0    5.0
 1.0  4.0   9.0  16.0   25.0
 1.0  8.0  27.0  64.0  125.0
```

Alternatively it is also possible to specify `obsdim = 1`, which
will cause each row to denote an observation.

```@jldoctest
julia> expand_poly(1:5, degree = 3, obsdim = 1)
5×3 Array{Int64,2}:
 1   1    1
 2   4    8
 3   9   27
 4  16   64
 5  25  125

julia> expand_poly(1:5, 3, ObsDim.First()); # same but type-stable
```
"""
function expand_poly(x::AbstractVector; degree::Integer = 5, obsdim = ObsDim.Last())
    expand_poly(x, degree, convert(LearnBase.ObsDimension, obsdim))
end

function expand_poly(x::AbstractVector, degree::Integer, ::ObsDim.Last)
    expand_poly(x, degree, ObsDim.Constant{2}())
end

function expand_poly(
    x::AbstractVector{T},
    degree::Integer,
    ::ObsDim.Constant{2} = ObsDim.Constant{2}()
) where T

    n = length(x)
    X = zeros(floattype(T), (degree, n))
    for i in 1:n
        for d in 1:degree
            @inbounds X[d, i] += float(x[i])^d
        end
    end
    return X
end

function expand_poly(
    x::AbstractVector{T},
    degree::Integer,
    ::ObsDim.Constant{1}
) where T

    n = length(x)
    X = zeros(floattype(T), (n, degree))
    for d in 1:degree
        for i in 1:n
            @inbounds X[i, d] += float(x[i])^d
        end
    end
    return X
end

floattype(::Type{T}) where T = float(T)
floattype(::Type{Missing}) = Missing
floattype(U::Union) = Union{floattype(U.a), floattype(U.b)}


text:
using Test
using Dates

include("gigasecond.jl")

samples = Dict(
    DateTime("2011-04-25") => DateTime("2043-01-01T01:46:40"),
    DateTime("1977-06-13") => DateTime("2009-02-19T01:46:40"),
    DateTime("1959-07-19") => DateTime("1991-03-27T01:46:40"),
    DateTime("2015-01-24T22:00:00") => DateTime("2046-10-02T23:46:40"),
    DateTime("2015-01-24T23:59:59") => DateTime("2046-10-03T01:46:39")
)

@testset "add gigasecond to $sample[1]" for sample in samples
    @test add_gigasecond(sample[1]) == sample[2]
end


text:
"""
    validate(mp::ONNX.ModelProto)
    validate(mp::ONNX.ModelProto, fs...)

Validate `mp`, throwing an exception if it is invalid.

It is possible to specify the validation steps `fs` to perform. Default is `uniqueoutput, optypedefined, outputused, inputused, hasname`
"""
validate(mp::ONNX.ModelProto, fs...=(uniqueoutput, optypedefined, outputused, inputused, hasname)...) = foreach(f -> f(mp), fs)

errinfo(n::ONNX.NodeProto) = "NodeProto with: \n"  * join(["\t$pn:\t$(clstring(getproperty(n, pn)))" for pn in propertynames(n) if hasproperty(n, pn)], "\n")
clstring(x::AbstractArray) = "[" * join(string.(x), ", ") * "]"
clstring(x) = string(x)

"""
    uniqueoutput(mp::ONNX.ModelProto, or=error)
    uniqueoutput(gp::ONNX.GraphProto, or=error)

Test that output names are unique. If not, an error message will be passed to `or`.
"""
uniqueoutput(mp::ONNX.ModelProto, or=error) = uniqueoutput(mp.graph, or)
function uniqueoutput(gp::ONNX.GraphProto, or=error)
    d = Dict()
    for n in gp.node
        for oname in n.output
            if haskey(d, oname)
                or("Duplicate output name: $oname found in \n$(errinfo(d[oname])) \nand\n $(errinfo(n))")
            end
            d[oname] = n
        end
    end
end

"""
    optypedefined(mp::ONNX.ModelProto, or=error)
    optypedefined(gp::ONNX.GraphProto, or=error)

Test that operations are defined for each node. If not, an error message will be passed to `or`.
"""
optypedefined(mp::ONNX.ModelProto, or=error) = optypedefined(mp.graph, or)
function optypedefined(gp::ONNX.GraphProto, or=error)
    for n in gp.node
        hasproperty(n, :op_type) || or("No op_type defined for $(errinfo(n))")
    end
end

"""
    outputused(mp::ONNX.ModelProto, or=error)
    outputused(gp::ONNX.GraphProto, or=error)

Test that all outputs are used. If not, an error message will be passed to `or`.
"""
outputused(mp::ONNX.ModelProto, or=error) = outputused(mp.graph, or)
function outputused(gp::ONNX.GraphProto, or=error)
    found, used = ioused(gp)
    unusedouts = setdiff(found, used)
    str(s) = join(sort(collect(s)), ", ")
    isempty(unusedouts) || or("Found unused outputs: $(str(unusedouts))")
end

"""
    inputused(mp::ONNX.ModelProto, or=error)
    inputused(gp::ONNX.GraphProto, or=error)

Test that all inputs are used. If not, an error message will be passed to `or`.
"""
inputused(mp::ONNX.ModelProto, or=error) = inputused(mp.graph, or)
function inputused(gp::ONNX.GraphProto, or=error)
    used, found = ioused(gp)
    unusedins = setdiff(found, used)
    str(s) = join(sort(collect(s)), ", ")
    isempty(unusedins) || or("Found unused inputs: $(str(unusedins))")
end

function ioused(gp::ONNX.GraphProto)
    found = union(Set(name.(gp.input)), Set(name.(gp.initializer)))
    used = Set(name.(gp.output))
    for n in gp.node
        foreach(oname -> push!(found, oname), n.output)
        foreach(iname -> push!(used, iname), n.input)
    end

    return found, used
end

hasname(mp::ONNX.ModelProto, or=error) = hasname(mp.graph, or)
function hasname(gp::ONNX.GraphProto, or=error)
     hasproperty(gp, :name) || return or("Graph name not defined!")
     isempty(gp.name) && or("Graph name is empty string!")
 end


text:
function get_projection(obj::BoundaryProjector)
    return jcall(obj, "getProjection", BoundaryProjection, ())
end

function visit_internal_node(obj::BoundaryProjector, arg0::BSPTree)
    return jcall(obj, "visitInternalNode", void, (BSPTree,), arg0)
end

function visit_leaf_node(obj::BoundaryProjector, arg0::BSPTree)
    return jcall(obj, "visitLeafNode", void, (BSPTree,), arg0)
end

function visit_order(obj::BoundaryProjector, arg0::BSPTree)
    return jcall(obj, "visitOrder", BSPTreeVisitor_Order, (BSPTree,), arg0)
end


text:
include("test_base.jl")

@testset "MLFlow" begin
    mlf = MLFlow()
    @test mlf.baseuri == "http://localhost:5000"
    @test mlf.apiversion == 2.0
    mlf = MLFlow("https://localhost:5001", apiversion=3.0)
    @test mlf.baseuri == "https://localhost:5001"
    @test mlf.apiversion == 3.0
end

@testset "createexperiment" begin
    @ensuremlf
    exp = createexperiment(mlf)
    @test isa(exp, MLFlowExperiment)
    @test deleteexperiment(mlf, exp)
    experiment = getexperiment(mlf, exp.experiment_id)
    @test experiment.experiment_id == exp.experiment_id
    @test experiment.lifecycle_stage == "deleted"
end

@testset "createrun" begin
    @ensuremlf
    expname = "getorcreate-$(UUIDs.uuid4())"
    e = getorcreateexperiment(mlf, expname)
    r = createrun(mlf, e.experiment_id)
    @test isa(r, MLFlowRun)
    @test deleterun(mlf, r)
    rr = createrun(mlf, e)
    @test isa(rr, MLFlowRun)
    @test deleterun(mlf, rr)
    @test deleteexperiment(mlf, e)
end

@testset "getorcreateexperiment" begin
    @ensuremlf
    expname = "getorcreate"
    e = getorcreateexperiment(mlf, expname)
    @test isa(e, MLFlowExperiment)
    ee = getorcreateexperiment(mlf, expname)
    @test isa(ee, MLFlowExperiment)
    @test e === ee
    @test deleteexperiment(mlf, ee)
    @test deleteexperiment(mlf, ee)
end

@testset "generatefilterfromparama" begin
    filter_params = Dict("k1" => "v1")
    filter = generatefilterfromparams(filter_params)
    @test filter == "param.\"k1\" = \"v1\""
    filter_params = Dict("k1" => "v1", "started" => Date("2020-01-01"))
    filter = generatefilterfromparams(filter_params)
    @test occursin("param.\"k1\" = \"v1\"", filter)
    @test occursin("param.\"started\" = \"2020-01-01\"", filter)
    @test occursin(" and ", filter)
end

@testset "searchruns" begin
    @ensuremlf
    exp = createexperiment(mlf)
    expid = exp.experiment_id
    exprun = createrun(mlf, exp)
    @test exprun.info.experiment_id == expid
    @test exprun.info.lifecycle_stage == "active"
    @test exprun.info.status == MLFlowRunStatus("RUNNING")
    exprunid = exprun.info.run_id

    runparams = Dict(
        "k1" => "v1",
        "started" => Date("2020-01-01")
    )
    logparam(mlf, exprun, runparams)

    findrun = searchruns(mlf, exp; filter_params=runparams)
    @test length(findrun) == 1
    r = only(findrun)
    @test get_run_id(get_info(r)) == exprun.info.run_id
    @test get_run_id(r) == get_run_id(get_info(r))
    @test sort(collect(keys(get_params(get_data(r))))) == sort(string.(keys(runparams)))
    @test sort(collect(values(get_params(get_data(r))))) == sort(string.(values(runparams)))
    @test get_params(r) == get_params(get_data(r))
    @test deleteexperiment(mlf, exp)
end

@testset "artifacts" begin
    @ensuremlf
    exp = createexperiment(mlf)
    @test isa(exp, MLFlowExperiment)
    exprun = createrun(mlf, exp)
    @test isa(exprun, MLFlowRun)
    # only run the below if artifact_uri is a local directory
    # i.e. when running mlflow server as a separate process next to the testset
    # when running mlflow in a container, the below tests will be skipped
    # this is what happens in github actions - mlflow runs in a container, the artifact_uri is not immediately available, and tests are skipped
    artifact_uri = exprun.info.artifact_uri
    if isdir(artifact_uri)
        @test_throws SystemError logartifact(mlf, exprun, "/etc/shadow")

        tmpfiletoupload = "sometempfilename.txt"
        f = open(tmpfiletoupload, "w")
        write(f, "samplecontents")
        close(f)
        artifactpath = logartifact(mlf, exprun, tmpfiletoupload)
        @test isfile(artifactpath)
        rm(tmpfiletoupload)
        artifactpath = logartifact(mlf, exprun, "randbytes.bin", b"some rand bytes here")
        @test isfile(artifactpath)

        mkdir(joinpath(artifact_uri, "newdir"))
        artifactpath = logartifact(mlf, exprun, joinpath("newdir", "randbytesindir.bin"), b"bytes here")
        artifactpath = logartifact(mlf, exprun, joinpath("newdir", "randbytesindir2.bin"), b"bytes here")
        mkdir(joinpath(artifact_uri, "newdir", "new2"))
        artifactpath = logartifact(mlf, exprun, joinpath("newdir", "new2", "randbytesindir.bin"), b"bytes here")
        artifactpath = logartifact(mlf, exprun, joinpath("newdir", "new2", "randbytesindir2.bin"), b"bytes here")
        mkdir(joinpath(artifact_uri, "newdir", "new2", "new3"))
        artifactpath = logartifact(mlf, exprun, joinpath("newdir", "new2", "new3", "randbytesindir.bin"), b"bytes here")
        artifactpath = logartifact(mlf, exprun, joinpath("newdir", "new2", "new3", "randbytesindir2.bin"), b"bytes here")
        mkdir(joinpath(artifact_uri, "newdir", "new2", "new3", "new4"))
        artifactpath = logartifact(mlf, exprun, joinpath("newdir", "new2", "new3", "new4", "randbytesindir.bin"), b"bytes here")
        artifactpath = logartifact(mlf, exprun, joinpath("newdir", "new2", "new3", "new4", "randbytesindir2.bin"), b"bytes here")

        # artifact tree should now look like this:
        #
        # ├── newdir
        # │   ├── new2
        # │   │   ├── new3
        # │   │   │   ├── new4
        # │   │   │   │   ├── randbytesindir2.bin
        # │   │   │   │   └── randbytesindir.bin
        # │   │   │   ├── randbytesindir2.bin
        # │   │   │   └── randbytesindir.bin
        # │   │   ├── randbytesindir2.bin
        # │   │   └── randbytesindir.bin
        # │   ├── randbytesindir2.bin
        # │   └── randbytesindir.bin
        # ├── randbytes.bin
        # └── sometempfilename.txt

        # 4 directories, 10 files

        artifactlist = listartifacts(mlf, exprun)
        @test sort(basename.(get_path.(artifactlist))) == ["newdir", "randbytes.bin", "sometempfilename.txt"]
        @test sort(get_size.(artifactlist)) == [0, 14, 20]

        ald2 = listartifacts(mlf, exprun, maxdepth=2)
        @test length(ald2) == 6
        @test sort(basename.(get_path.(ald2))) == ["new2", "newdir", "randbytes.bin", "randbytesindir.bin", "randbytesindir2.bin", "sometempfilename.txt"]
        aldrecursion = listartifacts(mlf, exprun, maxdepth=-1)
        @test length(aldrecursion) == 14 # 4 directories, 10 files
        @test sum(typeof.(aldrecursion) .== MLFlowArtifactDirInfo) == 4 # 4 directories
        @test sum(typeof.(aldrecursion) .== MLFlowArtifactFileInfo) == 10 # 10 files
    end
    deleterun(mlf, exprun)
    deleteexperiment(mlf, exp)
end

@testset "MLFlowClient.jl" begin
    @ensuremlf
    exp = createexperiment(mlf)
    @test isa(exp, MLFlowExperiment)

    exptags = [:key => "val"]
    expname = "expname-$(UUIDs.uuid4())"

    @test ismissing(getexperiment(mlf, "$(UUIDs.uuid4()) - $(UUIDs.uuid4())"))

    experiment = createexperiment(mlf; name=expname, tags=exptags)
    experiment_id = experiment.experiment_id
    experimentbyname = getexperiment(mlf, expname)
    @test experimentbyname.name == experiment.name

    exprun = createrun(mlf, experiment_id)
    @test exprun.info.experiment_id == experiment_id
    @test exprun.info.lifecycle_stage == "active"
    @test exprun.info.status == MLFlowRunStatus("RUNNING")
    exprunid = exprun.info.run_id

    logparam(mlf, exprunid, "paramkey", "paramval")
    logparam(mlf, exprunid, Dict("k" => "v", "k1" => "v1"))
    logparam(mlf, exprun, Dict("test1" => "test2"))

    logmetric(mlf, exprun, "metrickeyrun", 1.0)
    logmetric(mlf, exprun.info, "metrickeyrun", 2.0)
    logmetric(mlf, exprun.info, "metrickeyrun", [2.5, 3.5])
    logmetric(mlf, exprunid, "metrickey", 1.0)
    logmetric(mlf, exprunid, "metrickey2", [1.0, 1.5, 2.0])

    retrieved_run = getrun(mlf, exprunid)
    @test exprun.info == retrieved_run.info

    running_run = updaterun(mlf, exprunid, "RUNNING")
    @test running_run.info.experiment_id == experiment_id
    @test running_run.info.status == MLFlowRunStatus("RUNNING")
    finished_run = updaterun(mlf, exprun, MLFlowRunStatus("FINISHED"))
    finishedrun = getrun(mlf, finished_run.info.run_id)

    @test !ismissing(finishedrun.info.end_time)

    exprun2 = createrun(mlf, experiment_id)
    exprun2id = exprun.info.run_id
    logparam(mlf, exprun2, "param2", "key2")
    logmetric(mlf, exprun2, "metric2", [1.0, 2.0])
    updaterun(mlf, exprun2, "FINISHED")

    runs = searchruns(mlf, experiment_id)
    @test length(runs) == 2
    runs = searchruns(mlf, experiment_id; filter="param.param2 = \"key2\"")
    @test length(runs) == 1
    @test_throws ErrorException searchruns(mlf, experiment_id; run_view_type="ERR")
    runs = searchruns(mlf, experiment_id; filter="param.param2 = \"key3\"")
    @test length(runs) == 0
    runs = searchruns(mlf, experiment_id; max_results=1) # test paging functionality
    @test length(runs) == 2
    deleterun(mlf, exprunid)
    deleterun(mlf, exprun2)

    deleteexperiment(mlf, exp)
end

max_stars_repo_name: JuliaAI/MLFlowClient

max_issues_repo_name: JuliaAI/MLFlowClient

max_forks_repo_name: JuliaAI/MLFlowClient


text:
export minFrobNormspg

function minFrobNormspg(L::Tx, b, p, eta, indi, Lprev, weight, iter, fac)where {ETx<:Number, Tx<:AbstractMatrix{ETx}}
    ###### add for weighted part
    if weight!=1.0
        L*= weight
        L += (1/weight- 1) * Lprev * ( Lprev'*L)
    end

    ###### Matrix Dimension Setup
    undefine_1,m = size(b);
    undefine_2,r = size(L);
    R1 = Array{Any}(undef,m,r);

    ###### Solve for each row of R independently
        R1 = @distributed vcat for k = 1:length(indi)

            minnorm(L,b,k,eta, r, iter);

        end

        R = conj(permutedims(reshape(R1,r,m),[2,1]));

    return R
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule lttng_tools_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("lttng_tools")
JLLWrappers.@generate_main_file("lttng_tools", UUID("37591b3e-a3c2-50e0-b676-8bc36018f336"))
end  # module lttng_tools_jll


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Pod anti affinity is a group of inter pod anti affinity scheduling rules.

    IoK8sApiCoreV1PodAntiAffinity(;
        preferredDuringSchedulingIgnoredDuringExecution=nothing,
        requiredDuringSchedulingIgnoredDuringExecution=nothing,
    )

    - preferredDuringSchedulingIgnoredDuringExecution::Vector{IoK8sApiCoreV1WeightedPodAffinityTerm} : The scheduler will prefer to schedule pods to nodes that satisfy the anti-affinity expressions specified by this field, but it may choose a node that violates one or more of the expressions. The node that is most preferred is the one with the greatest sum of weights, i.e. for each node that meets all of the scheduling requirements (resource request, requiredDuringScheduling anti-affinity expressions, etc.), compute a sum by iterating through the elements of this field and adding \"weight\" to the sum if the node has pods which matches the corresponding podAffinityTerm; the node(s) with the highest sum are the most preferred.
    - requiredDuringSchedulingIgnoredDuringExecution::Vector{IoK8sApiCoreV1PodAffinityTerm} : If the anti-affinity requirements specified by this field are not met at scheduling time, the pod will not be scheduled onto the node. If the anti-affinity requirements specified by this field cease to be met at some point during pod execution (e.g. due to a pod label update), the system may or may not try to eventually evict the pod from its node. When there are multiple elements, the lists of nodes corresponding to each podAffinityTerm are intersected, i.e. all terms must be satisfied.
"""
mutable struct IoK8sApiCoreV1PodAntiAffinity <: SwaggerModel
    preferredDuringSchedulingIgnoredDuringExecution::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1WeightedPodAffinityTerm} } # spec name: preferredDuringSchedulingIgnoredDuringExecution
    requiredDuringSchedulingIgnoredDuringExecution::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PodAffinityTerm} } # spec name: requiredDuringSchedulingIgnoredDuringExecution

    function IoK8sApiCoreV1PodAntiAffinity(;preferredDuringSchedulingIgnoredDuringExecution=nothing, requiredDuringSchedulingIgnoredDuringExecution=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodAntiAffinity, Symbol("preferredDuringSchedulingIgnoredDuringExecution"), preferredDuringSchedulingIgnoredDuringExecution)
        setfield!(o, Symbol("preferredDuringSchedulingIgnoredDuringExecution"), preferredDuringSchedulingIgnoredDuringExecution)
        validate_property(IoK8sApiCoreV1PodAntiAffinity, Symbol("requiredDuringSchedulingIgnoredDuringExecution"), requiredDuringSchedulingIgnoredDuringExecution)
        setfield!(o, Symbol("requiredDuringSchedulingIgnoredDuringExecution"), requiredDuringSchedulingIgnoredDuringExecution)
        o
    end
end # type IoK8sApiCoreV1PodAntiAffinity

const _property_map_IoK8sApiCoreV1PodAntiAffinity = Dict{Symbol,Symbol}(Symbol("preferredDuringSchedulingIgnoredDuringExecution")=>Symbol("preferredDuringSchedulingIgnoredDuringExecution"), Symbol("requiredDuringSchedulingIgnoredDuringExecution")=>Symbol("requiredDuringSchedulingIgnoredDuringExecution"))
const _property_types_IoK8sApiCoreV1PodAntiAffinity = Dict{Symbol,String}(Symbol("preferredDuringSchedulingIgnoredDuringExecution")=>"Vector{IoK8sApiCoreV1WeightedPodAffinityTerm}", Symbol("requiredDuringSchedulingIgnoredDuringExecution")=>"Vector{IoK8sApiCoreV1PodAffinityTerm}")
Base.propertynames(::Type{ IoK8sApiCoreV1PodAntiAffinity }) = collect(keys(_property_map_IoK8sApiCoreV1PodAntiAffinity))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodAntiAffinity }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PodAntiAffinity[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodAntiAffinity }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodAntiAffinity[property_name]

function check_required(o::IoK8sApiCoreV1PodAntiAffinity)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PodAntiAffinity }, name::Symbol, val)
end


text:
module Example1

export Point
struct Point{T}
    x::T
    y::T
end

function Base. +(p::Point, q::Point)
    Point(p.x+ q.x, p.y + q.y)
end

export MPoint
mutable struct MPoint{T}
    x::T
    y::T
end

greet() = println("Hello World!")

################################################################

export Polynomial
struct Polynomial{T}
    coeffs::Vector{T}
end

export evaluate
function evaluate(p::Polynomial,x)
    sum(c * x^(i-1) for (i,c) in enumerate(p.coeffs))
end

function Base. +(p::Polynomial, q::Polynomial)
    l = max(length(p.coeffs), length(q.coeffs))
    rcoeffs = zeros(l)
    for (i,c) in enumerate(p.coeffs)
        rcoeffs[i] +=c
    end
    Polynomial(rcoeffs)
end

function Base. *(a, p::Polynomial)
    Polynomial(map(x->a*x, p.coeffs))
end

#Base .*(p::Polynomial,a) = a*p

export deriv
function deriv(p::Polynomial)
    Polynomial([(i-1) * c for (i,c) in enumerate(p.coeffs)][2::end])
end

end # module


text:
using Benchmarks
using FixedSizeArrays

import Base: *, isless, isequal

include("utils.jl");
include("algorithms.jl");

function graham_scan_vanilla(p::Vector{Point{2, Float64}})
    ep = _grahamscan(p)

    return ep
end

function graham_scan_wat(p::Vector{Point{2, Float64}})
    p = _akltoussaint(p)
    ep = _grahamscan(p)

    return ep
end

function monotonechain_vanilla(p::Vector{Point{2, Float64}})
    ep = _monotonechain(p)

    return ep
end

function monotonechain_wat(p::Vector{Point{2, Float64}})
    p = _akltoussaint(p)
    ep = _monotonechain(p)

    return ep
end

p = [Point(randn(), randn()) for i=1:2500000]

srand(1112016)
npoints = [25, 250, 2500, 25000, 250000, 2500000]
nreps = [1000, 500, 100, 50, 25, 5]

for i=1:length(npoints)

    # Get parameters
    np = npoints[i]
    nr = nreps[i]

    p = [Point(randn(), randn()) for j=1:np]

    gsv() = graham_scan_vanilla(p)
    gsat() = graham_scan_wat(p)
    mcv() = monotonechain_vanilla(p)
    mcat() = monotonechain_wat(p)

    println("With $np points")
    println(compare([gsv, gsat, mcv, mcat], nr))
end


text:
"""
function plotNodeVisited(metrics::BasicMetrics{O, H}; filename::String="") where{O<:AbstractTakeObjective, H<:ValueSelection}

plot 2 graphs :
1) number of node visited by the heuristic to find a first solution for every learning episode.
2) number of node visited by the heuristic to prove to optimality for every learning episode.

The learning process should show a decrease in the number of nodes required to find a first solution along the search
(depending on the reward engineering).
"""
function plotNodeVisited(metricsArray::Union{BasicMetrics, Vector{AbstractMetrics}}; filename::String="")
    L = length(metricsArray[1].meanNodeVisitedUntilEnd)
    Label = Matrix{String}(undef, 1, length(metricsArray))
    learnedIdx = 1
    basicIdx = 1
    for j in 1:length(metricsArray) #nb of heuristics
        if isa(metricsArray[j].heuristic, LearnedHeuristic)
            Label[1,j]="Learned heuristic n°$learnedIdx"
            learnedIdx += 1
        else
            Label[1,j]="Classic heuristic n°$basicIdx"
            basicIdx += 1
        end
    end
    p1 = plot(
        1:L,
        [metricsArray[i].meanNodeVisitedUntilEnd[1:L] for i in 1:length(metricsArray)],
        title = "Node visited until optimality",
        label=Label,
        xlabel="Episode",
        ylabel="Nodes visited",
    )
    p2 = plot(
        1:L,
        [metricsArray[i].meanNodeVisitedUntilfirstSolFound[1:L] for i in 1:length(metricsArray)],
        title ="Node visited until first solution found",
        label=Label,
        xlabel="Episode",
        ylabel="Nodes visited",
    )
    plot(p1,p2, layout=(2,1))
end

function plotNodeVisited(metrics::BasicMetrics; filename::String="")
    L = length(metrics.meanNodeVisitedUntilEnd)
    p = plot(
        1:L,
        [metrics.meanNodeVisitedUntilEnd[1:L] metrics.meanNodeVisitedUntilfirstSolFound[1:L]],
        xlabel="Episode",
        ylabel="Nodes visited",
        title = ["Node visited until Optimality" "Node visited until first solution found"],
        layout = (2, 1)
    )
    display(p)
    savefig(p,filename*"_node_visited_"*"$(typeof(metrics.heuristic))"*".png")
end

"""
function plotScoreVariation(metrics::BasicMetrics{O, H}; filename::String="") where{O<:AbstractTakeObjective, H<:ValueSelection}

plot the relative scores ( compared to the optimal ) of the heuristic during the search for fixed instances along the training. This plot is
meaningful only if the metrics is one from the evaluator (ie. the instance remains the same one).
"""
function plotScoreVariation(metrics::BasicMetrics{TakeObjective, H}; filename::String="") where H
    Data=[]
    for i in length(metrics.nodeVisited):-1:1
        push!(Data,hcat(metrics.nodeVisited[i],metrics.scores[i]))
    end

    p = plot(
        [scatter[:,1] for scatter in Data],
        [scatter[:,2] for scatter in Data],
        #fillrange =[[ones(length(scatter[:,2])),scatter[:,2]] for scatter in Data],
        #fillalpha=1,
        xlabel="number of nodes visited",
        ylabel="relative score",
        xaxis=:log,
    )
    display(p)
    savefig(p,filename*"_score_variation.png")
end

"""
function plotScoreVariation(metrics::BasicMetrics{O, H}; filename::String="") where{O<:AbstractTakeObjective, H<:ValueSelection}

plot the relative scores ( compared to the optimal ) of the heuristic during the search for fixed instances along the training. This plot is
meaningful only if the metrics is one from the evaluator (ie. the instance remains the same one).
"""
function plotRewardVariation(metrics::BasicMetrics{<:AbstractTakeObjective, <:LearnedHeuristic}; filename::String="")
    L = length(metrics.totalReward)

    p = plot(1:L,
        metrics.totalReward[1:L],
        xlabel="Episode",
        ylabel="Total reward",
        title ="Total reward per episode",
    )
    display(p)
    savefig(p,filename*"_reward_variation.png")
end


text:
@time @testset "Properties.Property" begin include("unit_Property.jl") end
@time @testset "Properties.tune" begin include("unit_tune.jl") end


text:
number_of_mRNA_species_balances = 23
number_of_protein_species_balances = 23
number_of_auxilary_balances = 0
gene_deletion_array = ones(number_of_mRNA_species_balances)

# Assume we have no mRNA or protein initially, but we do have cells -
total_number_of_species = number_of_mRNA_species_balances+number_of_mRNA_species_balances+number_of_auxilary_balances
initial_condition_vector = zeros(total_number_of_species)
cell_mass_index = number_of_mRNA_species_balances+number_of_mRNA_species_balances+1
atra_index = number_of_mRNA_species_balances+number_of_mRNA_species_balances+2
cRAFS621_index = number_of_mRNA_species_balances+number_of_mRNA_species_balances+3

# constants (from bionumbers)
hl60_diameter = 12                # \mum
number_of_rnapII = 75000          # copies/cells
number_of_ribosome = 1e6          # copies/cells
mRNA_half_life_TF = 2             # hrs
protein_half_life = 10            # hrs
doubling_time_hl60 = 19.5         # hrs
max_translation_rate = 5          # aa/sec
max_transcription_rate = 6.0       # n/sec
average_transcript_length = 15000 # n
average_protein_length = 5000     # aa
fraction_nucleus = 0.49           # dimensionless
av_number = 6.02e23               # number/mol
avg_gene_number = 2               # number of copues of a gene -

# Calculate the volume (convert to L)
V = ((1-fraction_nucleus)*(1/6)*(3.14159)*(hl60_diameter)^3)*(1e-15)

# Calculate the rnapII_concentration and ribosome_concentration
rnapII_concentration = number_of_rnapII*(1/av_number)*(1/V)*1e9 # \nM
ribosome_concentration = number_of_ribosome*(1/av_number)*(1/V)*1e9 # \nM

# degrdation rate constants -
degradation_constant_mRNA = -(1/mRNA_half_life_TF)*log(0.5)
degradation_constant_protein = -(1/protein_half_life)*log(0.5)

# kcats for transcription and translation -
kcat_transcription = max_transcription_rate*(3600/average_transcript_length)  # hr^-1
kcat_translation = max_translation_rate*(3600/average_protein_length)  # hr^-1

# Maximum specific growth rate -
maximum_specific_growth_rate = (1/doubling_time_hl60)*log(2) # 0.0355 # hr^-1

# What is the average gene concentration -
avg_gene_concentration = avg_gene_number*(1/av_number)*(1/V)*1e9  #\nM

# What fraction of the maximum is the background?
background_expression_parameter_array = 0.0000001*ones(number_of_mRNA_species_balances)
background_expression_parameter_array[1] = 0.0
background_expression_parameter_array[2] = 0.0
background_expression_parameter_array[3] = 0.0

# How fast do my ells die?
death_rate_constant = 0.2*maximum_specific_growth_rate

# Saturation constants for translation and trascription -
saturation_transcription = 4600*(1/av_number)*(1/V)*1e9
saturation_translation = 100000*(1/av_number)*(1/V)*1e9

const_transcription = kcat_transcription/saturation_transcription
const_translation   = kcat_translation/saturation_translation

eg_transciption_rate = 1.0 * const_transcription * rnapII_concentration
eg_translation_rate = 1.0 * const_translation * ribosome_concentration


text:
using Mimi, CSVFiles, DataFrames, Query, Interpolations

@defcomp SSPs begin

    country = Index()

    SSP_source          = Parameter{String}() # can be one of IIASA GDP, OECD Env-Growth, PIK GDP_32, and Benveniste
    SSP                 = Parameter{String}() # can be one of SSP1, SSP2, SSP3, SSP4, SSP5
    emissions_source    = Parameter{String}() # can be one of Leach
    emissions_scenario  = Parameter{String}() # can be one of SSP119, SSP126, SSP245, SSP370, SSP585

    country_names       = Parameter{String}(index=[country]) # need the names of the countries from the dimension

    population          = Variable(index=[time, country], unit="million")
    population_global   = Variable(index=[time], unit="million")
    gdp                 = Variable(index=[time, country], unit="billion US\$2005/yr")
    gdp_global          = Variable(index=[time], unit="billion US\$2005/yr")

    co2_emissions       = Variable(index=[time], unit="GtC/yr")
    ch4_emissions       = Variable(index=[time], unit="MtCH4/yr")
    n2o_emissions       = Variable(index=[time], unit="MtN/yr")
    sf6_emissions       = Variable(index=[time], unit="MtSF6/yr")

    function init(p,v,d)

        # ----------------------------------------------------------------------
        # Checks

        ssp_model_options = ["IIASA GDP", "OECD Env-Growth", "PIK GDP_23", "Benveniste"]
        !(p.SSP_source in ssp_model_options) && error("Model $(p.SSP_source) provided to SSPs component SSP_source parameter not found in available list: $(ssp_model_options)")

        ssp_options = ["SSP1", "SSP2", "SSP3", "SSP4", "SSP5"]
        !(p.SSP in ssp_options) && error("SSP $(p.SSP) provided to SSPs component SSP parameter not found in available list: $(ssp_options)")

        emissions_source_options = ["Leach"]
        !(p.emissions_source in emissions_source_options) && error("Model $(p.emissions_source) provided to SSPs component emissions_source parameter not found in available list: $(emissions_source_options)")

        emissions_scenario_options = ["SSP119", "SSP126", "SSP245", "SSP370", "SSP585"]
        !(p.emissions_scenario in emissions_scenario_options) && error("emissions_scenario $(p.emissions_scenario) provided to SSPs component emissions_scenario parameter not found in available list: $(emissions_scenario_options)")

        # ----------------------------------------------------------------------
        # Load Socioeconomic Data as Needed
        #   population in millions of individuals
        #   GDP in billions of $2005 USD

        socioeconomic_path = joinpath(@__DIR__, "..", "..", "data", "socioeconomic", "$(p.SSP_source)_$(p.SSP).csv")
        ssp_dict_key = Symbol(p.SSP_source, "-", p.SSP)

        if !haskey(g_ssp_datasets, ssp_dict_key)
            g_ssp_datasets[ssp_dict_key] = load(socioeconomic_path) |> DataFrame
        end

        # Check Countries - each country found in the model countries parameter
        # must exist in the SSP socioeconomics dataframe

        missing_countries = []
        unique_country_list = unique(g_ssp_datasets[ssp_dict_key].country)
        for country in p.country_names
            !(country in unique_country_list) && push!(missing_countries, country)
        end
        !isempty(missing_countries) && error("All countries in countries parameter must be found in SSPs component Socioeconomic Dataframe, the following were not found: $(missing_countries)")

        # ----------------------------------------------------------------------
        # Load Emissions Data as Needed
        #   carbon dioxide emissions in GtC
        #   nitrous oxide emissions in MtN
        #   methane emissions in MtCH4
        #   SF6 emissions in MtSF6

        emissions_path = joinpath(@__DIR__, "..", "..", "data", "emissions", "$(p.emissions_source)_$(p.emissions_scenario).csv")
        emissions_scenario_dict_key = Symbol(p.emissions_source, "-", p.emissions_scenario)

        if !haskey(g_emissions_scenario_datasets, emissions_scenario_dict_key)
            if p.emissions_source == "Leach"
                emissions_data = load(emissions_path, skiplines_begin = 6)|>
                        DataFrame |>
                        i -> rename!(i, Symbol.([:year, names(i)[2:end]...])) |>
                        DataFrame |>
                        @select(:year, :carbon_dioxide, :nitrous_oxide, :methane, :sf6) |>
                        DataFrame
            else # already perfect formatted
                emissions_data = load(emissions_path)|> DataFrame
            end

            g_emissions_scenario_datasets[emissions_scenario_dict_key]= emissions_data
        end
    end

    function run_timestep(p,v,d,t)

        ssp_dict_key = Symbol(p.SSP_source, "-", p.SSP)
        emissions_scenario_dict_key = Symbol(p.emissions_source, "-", p.emissions_scenario)

        year_label = gettime(t)

        # check that we only run the component where we have data
        if !(year_label in g_ssp_datasets[ssp_dict_key].year)
            error("Cannot run SSP component in year $(year_label), SSP socioeconomic variables not available for this model and year.")
        end
        if !(year_label in g_emissions_scenario_datasets[emissions_scenario_dict_key].year)
            error("Cannot run SSP component in year $(year_label), SSP emissions variables not available for this model and year.")
        end

        # ----------------------------------------------------------------------
        # Socioeconomic

        # filter the dataframe for values with the year matching timestep
        # t and only the SSP countries found in the model countries list,
        # already checked that all model countries are in SSP countries list
        subset = g_ssp_datasets[ssp_dict_key] |>
            @filter(_.year == gettime(t) && _.country in p.country_names) |>
            DataFrame

        # get the ordered indices of the SSP countries within the parameter
        # of the model countries, already checked that all model countries
        # are in SSP countries list
        order = indexin(p.country_names, subset.country)

        v.population[t,:] = subset.pop[order]
        v.gdp[t,:] = subset.gdp[order]

        # add global data for future accessibility and quality control
        v.population_global[t] = sum(v.population[t,:])
        v.gdp_global[t] = sum(v.gdp[t,:])

        # ----------------------------------------------------------------------
        # Emissions

        subset = g_emissions_scenario_datasets[emissions_scenario_dict_key] |>
                    @filter(_.year == gettime(t)) |>
                    DataFrame

        v.co2_emissions[t] = subset.carbon_dioxide[1]
        v.ch4_emissions[t] = subset.methane[1]
        v.n2o_emissions[t] = subset.nitrous_oxide[1]
        v.sf6_emissions[t] = subset.sf6[1]

    end
end


text:
using Makie
using MakieLayout
using Animations

begin
    container_scene = Scene(camera = campixel!, resolution = (1200, 1200))
    # display(container_scene)

    t = Node(0.0)

    a_width = Animation([1, 7], [1200, 800], sineio(n=2, yoyo=true, postwait=0.5))
    a_height = Animation([2.5, 8.5], [1200, 800], sineio(n=2, yoyo=true, postwait=0.5))

    scene_area = lift(t) do t
        IRect(0, 0, round(Int, a_width(t)), round(Int, a_height(t)))
    end

    scene = Scene(container_scene, scene_area, camera = campixel!)

    rect = poly!(scene, scene_area,
        raw=true, color=RGBf0(0.97, 0.97, 0.97), strokecolor=:transparent, strokewidth=0)[end]

    outer_gl = GridLayout(scene, alignmode = Outside(30))

    inner_gl = outer_gl[1, 1] = GridLayout()

    ax1 = inner_gl[1, 1] = LAxis(scene, xautolimitmargin=(0, 0), yautolimitmargin=(0, 0))
    ax2 = inner_gl[1, 2] = LAxis(scene, xautolimitmargin=(0, 0), yautolimitmargin=(0, 0))
    ax3 = inner_gl[2, 1:2] = LAxis(scene, xautolimitmargin=(0, 0), yautolimitmargin=(0, 0))

    guigl = inner_gl[3, 1:2] = GridLayout()
    b1 = guigl[1, 1] = LButton(scene, label = "prev", width = Auto())
    sl = guigl[1, 2] = LSlider(scene, startvalue = 6, height = 40)
    b2 = guigl[1, 3] = LButton(scene, label = "next", width = Auto())

    data = randn(200, 200) .+ 3 .* sin.((1:200) ./ 20) .* sin.((1:200)' ./ 20)
    h1 = heatmap!(ax1, data)
    h2 = heatmap!(ax2, data, colormap = :blues)
    h3 = heatmap!(ax3, data, colormap = :heat)

    agl1 = gridnest!(inner_gl, 1, 1)
    agl1[1, 2] = LColorbar(scene, h1, width = 30, label = "normal bar")
    agl1[2, 1:2] = LSlider(scene, height = 20, startvalue = 4)
    agl1[3, 1:2] = LSlider(scene, height = 20, startvalue = 5)
    agl1[4, 1:2] = LSlider(scene, height = 20, startvalue = 6)

    agl2 = gridnest!(inner_gl, 1, 2)
    agl2[1, 2] = LColorbar(scene, h2, width = 30, height = Relative(0.66), label = "two thirds bar")
    agl2gl = agl2[2, :] = GridLayout()
    agl2gl[1, 1] = LButton(scene, label = "Run", height = Auto())
    agl2gl[1, 2] = LButton(scene, label = "Start")

    agl3 = gridnest!(inner_gl, 2, 1:2)
    agl3[:, 3] = LColorbar(scene, h3, width = 30, height=200, label = "fixed height bar")
    rowsize!(agl3, 1, Auto(false, 1.0))

    inner_gl[0, :] = LText(scene, text = "MakieLayout", textsize = 50)
end

for ts in 0:1/30:10
    t[] = ts
    sleep(1/60)
end

record(container_scene, "/Users/juliuskrumbiegel/Desktop/layoutdemo.mp4", 0:1/30:9; framerate=30) do ti
    t[] = ti
end