text:
using LinearAlgebra
using Distributed
using SharedArrays
using REPL

export DSProblem, SetObjective, SetInitialPoint, SetVariableBound, SetMaxEvals,
       SetOpportunisticEvaluation, SetSense, SetVariableBounds, Optimize!, SetGranularity,
       SetGranularities

"""
	DSProblem{T}(N::Int; poll::AbstractPoll=LTMADS{T}(),
                         search::AbstractSearch=NullSearch(),
                         objective::Union{Function,Nothing}=nothing,
                         initial_point::Vector=zeros(T, N),
                         iteration_limit::Int=1000,
                         function_evaluation_limit::Int=5000,
                         sense::ProblemSense=Min,
                         full_output::Bool=false,
                         granularity::Vector=zeros(T, N),
                         min_mesh_size::Union{T,Nothing}=nothing,
                         min_poll_size::Union{T,Nothing}=nothing,
                         kwargs...
                         ) where T

Return a problem definition for an `N` dimensional problem.

`poll` and `search` specify the poll and search step algorithms to use. The default
choices are [`LTMADS`](@ref) and [`NullSearch`](@ref) respectively.

Note that if working with `Float64` (normally the case) then the type
parameterisation can be ignored.
"""
mutable struct DSProblem{T, MT, ST, PT, CT} <: AbstractProblem{T} where {MT <: AbstractMesh, ST <: AbstractSearch, PT <: AbstractPoll, CT <: AbstractCache}
    #= Problem Definition =#
    objective::Function
    constraints::Constraints{T}
    #Problem size
    N::Int
    user_initial_point::Union{Vector{T},Nothing}
    granularity::Vector{T}
    lower_bounds::Vector{Union{T, Nothing}}
    upper_bounds::Vector{Union{T, Nothing}}
    #Barrier threshold
    h_max::T
    sense::ProblemSense

    #TODO incumbent points should be sets not points, therefore change to vectors of points
    #and remove type unions

    #= Working Variables =#
    #Feasible incumbent point
    x::Union{Vector{T},Nothing}
    #Feasible incumbent point evaluated cost
    x_cost::Union{T,Nothing}
    #Infeasible incumbent point
    i::Union{Vector{T},Nothing}
    #Infeasible incumbent point evaluated cost
    i_cost::Union{T,Nothing}

    cache::CT

    #= Runtime data =#
    status::Status{T}

    #= Solver Config =#
    config::Config{T, MT, ST, PT}

    stoppingconditions::Vector{AbstractStoppingCondition}

    full_output::Bool

    DSProblem(N::Int;kwargs...) = DSProblem{Float64}(N; kwargs...)

    function DSProblem{T}(N::Int;
                          poll::AbstractPoll=LTMADS{T}(),
                          search::AbstractSearch=NullSearch(),
                          objective::Union{Function,Nothing}=nothing,
                          initial_point::Vector=zeros(T, N),
                          iteration_limit::Int=1000,
                          function_evaluation_limit::Int=5000,
                          sense::ProblemSense=Min,
                          full_output::Bool=false,
                          granularity::Vector=zeros(T, N),
                          min_mesh_size::Union{T,Nothing}=nothing,
                          min_poll_size::Union{T,Nothing}=nothing,
                          kwargs...
                         ) where T

        p = new{T, Mesh{T}, typeof(search), typeof(poll), PointCache{T}}()

        p.N = N
        p.user_initial_point = convert(Vector{T},initial_point)
        p.granularity = convert(Vector{T},granularity)
        p.lower_bounds = fill(nothing, N)
        p.upper_bounds = fill(nothing, N)

        p.sense = sense

        p.config = Config{T}(N, poll, search, Mesh{T}(N);kwargs...)

        p.status = Status{T}()
        p.cache = PointCache{T}()
        p.constraints = Constraints{T}()

        p.stoppingconditions = AbstractStoppingCondition[
            IterationStoppingCondition(iteration_limit),
            FunctionEvaluationStoppingCondition(function_evaluation_limit),
            MeshPrecisionStoppingCondition{T}(min_mesh_size),
            PollPrecisionStoppingCondition{T}(min_poll_size)
        ]

        p.x = nothing
        p.x_cost = nothing
        p.i = nothing
        p.i_cost = nothing

        if objective != nothing
            p.objective = objective
        end

        p.full_output = full_output

        return p
    end
end

MeshUpdate!(p::DSProblem, result::IterationOutcome) =
    MeshUpdate!(p.config.mesh, p.config.poll, result, p.status.success_direction)

"""
    SetMaxEvals(p::DSProblem, max::Bool=true)

Set/unset parallel blackbox evaluations.
The number of threads Julia was started with will be used.

Note that using parallel blackbox evaluations will only result in reduced runtime
for problems that have long blackbox evaluations.
"""
function SetMaxEvals(p::DSProblem, max::Bool=true)
    if max
        p.config.max_simultaneous_evaluations = p.config.num_threads
        if p.config.num_threads == 1
            println("Julia was started single-threaded.")
            println("Start Julia with the option `--threads N` where N is the number of threads,\n to use parallel blackbox evaluations.")
        end
    else
        p.config.max_simultaneous_evaluations = 1
    end
    p.config.parallel_lock = ReentrantLock()
end

"""
    SetObjective(p::DSProblem, obj::Function)

Sets the target objective function to solve. `obj` should take a vector and return
a single cost value.
"""
function SetObjective(p::DSProblem, obj::Function)
    if p.sense == Min
        p.objective = obj
    else
        p.objective = x -> -obj(x)
    end
end

"""
    SetSense(p::DSProblem, sense::ProblemSense )

Set the problem sense. Valid values for `sense` are `DS.Min` and `DS.Max`.
"""
function SetSense(p::DSProblem, sense::ProblemSense)
    p.sense = sense
end

"""
    SetInitialPoint(p::DSProblem{T}, x::Vector{T}) where T

Set the initial incumbent point to `x`. This must be of the correct dimension. If using
any extreme barrier constraints then it must also satisfy these constraints.
"""
function SetInitialPoint(p::DSProblem{T}, x::Vector{T}) where T
    size(x, 1) == p.N || error("Point dimensions don't match problem definition")
    #TODO Check against constraints
    p.user_initial_point = x
end

"""
    SetGranularity(p::DSProblem{T}, index::Int, g::T) where T

Set the granularity of the variable with index `i` to `g`.
"""
function SetGranularity(p::DSProblem{T}, index::Int, g::T) where T
    1 <= index <= p.N || error("Invalid variable index, should be in range 1 to $(p.N).")
    g >= 0 || error("Granularity has to be non-negative.")

    p.granularity[index] = g
end

"""
    SetGranularities(p::DSProblem{T}, g::Vector{T}) where T

Call [`SetGranularity`](@ref) for each variable. The vector `g` should contain the granularity
for each variable.
"""
function SetGranularities(p::DSProblem{T}, g::Vector{T}) where T
    size(g, 1) == p.N || error("Granularity vector dimensions don't match problem definition")

    for i=1:p.N
        SetGranularity(p, i, g[i])
    end
end

"""
    SetOpportunisticEvaluation(p::DSProblem; opportunistic::Bool=true)

Set/unset opportunistic evaluation (enables by default).

When using opportunistic evaluation the first allowable evaluated point with an
improved cost is set as the new incumbent solution. If using progressive barrier
constraints this point may be infeasible.
"""
function SetOpportunisticEvaluation(p::DSProblem; opportunistic::Bool=true)
    p.config.opportunistic = opportunistic
end

"""
    EvaluateInitialPoint(p::DSProblem)

Evaluate the initial point.

Throws an error if the initial point is not feasible.
"""
function EvaluateInitialPoint(p::DSProblem)
    p.user_initial_point == Nothing() && return
    feasibility = ConstraintEvaluation(p, p.user_initial_point)
    if feasibility == StrongInfeasible
        error("Initial point must be feasible")
    elseif feasibility == WeakInfeasible
        p.i = p.user_initial_point
        p.i_cost = round(p.objective(p.user_initial_point), digits=p.config.cost_digits)
        CachePush(p, p.i, p.i_cost)
    elseif feasibility == Feasible
        p.x = p.user_initial_point
        p.x_cost = round(p.objective(p.user_initial_point), digits=p.config.cost_digits)
        CachePush(p, p.x, p.x_cost)
    end

    p.status.function_evaluations += 1

    p.full_output && InitialPointEvaluationOutput(p, feasibility)
end

#TODO review if this is the kind of scaling we want to do
"""
	SetVariableBound(p::DSProblem{T}, index::Int, l::T, u::T) where T

Set the expected bounds of the variable with index `i` to between lower (`l`) and upper (`u`)
values. This **does not** create a constraint, and is only used for scaling when a variable
varies with a significantly different scale to the others.
"""
function SetVariableBound(p::DSProblem{T}, index::Int, l::T, u::T) where T
    1 <= index <= p.N || error("Invalid variable index, should be in range 1 to $(p.N).")

    p.lower_bounds[index] = isinf(l) ? nothing : l
    p.upper_bounds[index] = isinf(u) ? nothing : u
end

"""
	SetVariableBounds(p::DSProblem{T}, l::Vector{T}, u::Vector{T}) where T

Call [`SetVariableBound`](@ref) for each variable. The vectors `l` and `u` should contain
a lower and upper bound for each variable.
"""
function SetVariableBounds(p::DSProblem{T}, l::Vector{T}, u::Vector{T}) where T
    size(l, 1) == p.N || error("Lower bound vector dimensions don't match problem definition")
    size(u, 1) == p.N || error("Upper bound vector dimensions don't match problem definition")

    for i=1:p.N
        SetVariableBound(p, i, l[i], u[i])
    end
end

"""
    Optimize!(p::DSProblem)

Run the direct search algorithm on problem `p`.

`p` must have had its initial point and objective function set. If extreme
barrier constraints have been set then the initial point must be value for
those constraints.
"""

function Optimize!(p::DSProblem)
    # start to monitor the input in REPL (for KeyInterrupt)
    term = REPL.Terminals.TTYTerminal("xterm",stdin,stdout,stderr)
    REPL.Terminals.raw!(term,true)
    Base.start_reading(stdin)

    # check the dimension of the problem
    p_dimention=p_dim(p)
    if p_dimention>2
        error("Sorry, Multi-MADS are not available.")
    elseif p_dimention==2
        return Optimize_Bi!(p)
    else
    #TODO check that problem definition is complete
    Setup(p)

    while _check_stoppingconditions(p)
        p.full_output && OutputIterationDetails(p)
        OptimizeLoop(p)
    end

    Finish(p)
    end
end

#Initialise solver
function Setup(p)
    p.status.start_time = time()
    _check_initial_point(p)
    MeshSetup!(p)
    _init_stoppingconditions(p)
    EvaluateInitialPoint(p)
    CacheOrderPush(p)
end

#A single iteration of the search->poll->update algorithm loop
function OptimizeLoop(p)
    if p.full_output
        println("Search step:\n")
    end
    result = Search(p)

    if p.full_output && !(p.config.search isa NullSearch)
        println("\tResult of Search step: $result\n")
    end

    if p.full_output
        println("Poll step:\n")
        if result !== Unsuccessful
            println("\tSkipping Poll step.\n")
        end
    end

    #If search fails, run poll
    if result == Unsuccessful
        result = Poll(p)
        if p.full_output
            println("\tResult of Poll step: $result\n")
        end
    end

    result != Unsuccessful && CacheOrderPush(p)

    #pass the result of search/poll to update
    MeshUpdate!(p, result)

    p.status.iteration += 1
    return result
end

#Cleanup and reporting
function Finish(p)
    p.status.runtime_total = time() - p.status.start_time
    ReportFinal(p)
end

"""
    EvaluatePoint!(p::DSProblem{FT}, trial_points)::IterationOutcome where {FT<:AbstractFloat}

Determine whether the set of trial points result in a dominating, improving, or unsuccesful
algorithm iteration. Update the feasible and infeasible incumbent points of `p`.
"""
function EvaluatePoint!(p::DSProblem{FT}, trial_points::Vector{Vector{FT}})::IterationOutcome where {FT<:AbstractFloat}
    if p.config.max_simultaneous_evaluations > 1
        EvaluatePointParallel!(p, trial_points)
    else
        EvaluatePointSequential!(p, trial_points)
    end
end

"""
    valuatePointSequential!(p::DSProblem{FT}, trial_points::Vector{Vector{FT}})::IterationOutcome where {FT<:AbstractFloat}

Single-threaded evaluation of set of trial points.
"""
function EvaluatePointSequential!(p::DSProblem{FT}, trial_points::Vector{Vector{FT}})::IterationOutcome where {FT<:AbstractFloat}
    #TODO could split into an evaluation function and an update function
    isempty(trial_points) && return Unsuccessful

    #Variables to store the best evaluated points
    feasible_point = isnothing(p.x) ? FT(Inf) * ones(p.N) : p.x
    feasible_cost = isnothing(p.x_cost) ? FT(Inf) : p.x_cost
    infeasible_point = isnothing(p.i) ? FT(Inf) * ones(p.N) : p.i
    infeasible_cost = isnothing(p.i_cost) ? FT(Inf) : p.i_cost
    successful_direction = nothing

    #The current minimum hmax value for all collections
    h_min = GetOldHmaxSum(p.constraints)

    if p.full_output
        println("\tPoint evaluation:\n")
    end

    #Iterate over all trial points
    for (i, point)=enumerate(trial_points)

        feasibility = ConstraintEvaluation(p, point)

        # Point violates h_max on at least one constraint collection
        feasibility == StrongInfeasible && continue

        #Point is feasible for relaxed constraints, so evaluate it
        p.status.blackbox_time_total += @elapsed (cost, is_from_cache) = function_evaluation(p, point)

        #To determine if a point is dominating or improving the combined h_max is needed
        h = GetViolationSum(p.constraints, point)

        updated = false

        # Conditions met for a dominant point
        if feasibility == Feasible && cost < feasible_cost
            feasible_point = point
            feasible_cost = cost
            successful_direction = isnothing(p.status.directions) ? nothing : p.status.directions[i]
            updated = true
        elseif feasibility == WeakInfeasible && h < h_min
            # Conditions met for an improving point (worse cost, but closer to being feasible) or
            # a dominant point (better cost and closer to feasibility)
            # Only record if it offers an improved constraint violation
            infeasible_point = point
            infeasible_cost = cost
            h_min = h
            successful_direction = isnothing(p.status.directions) ? nothing : p.status.directions[i]
            updated = true
        end

        p.full_output && OutputPointEvaluation(i, point, cost, h, is_from_cache)

        #break if using opportunistic iteration
        updated && p.config.opportunistic && break
    end

    result = Unsuccessful

    incum_i_cost = isnothing(p.i_cost) ? FT(Inf) : p.i_cost
    incum_x_cost = isnothing(p.x_cost) ? FT(Inf) : p.x_cost

    p.status.directions = nothing

    # Dominates if there is a feasible improvement, or an infeasible point with
    # reduced violation (determined previously) as well as a cost lower than any
    # yet tested (feasible and infeasible)
    if feasible_cost < incum_x_cost
        result = Dominating
    elseif infeasible_cost < incum_i_cost && h_min < GetOldHmaxSum(p.constraints)
        result = Dominating
    elseif h_min < GetOldHmaxSum(p.constraints)
        result = Improving
    end

    if feasible_cost < incum_x_cost
        p.x = feasible_point
        p.x_cost = feasible_cost
    end

    if infeasible_cost < incum_i_cost || h_min < GetOldHmaxSum(p.constraints)
        p.i = infeasible_point
        p.i_cost = infeasible_cost
    end

    p.status.success_direction = successful_direction

    UpdateConstraints(p.constraints, h_min, result, p.x, p.i)

    return result
end

"""
    EvaluatePointParallel!(p::DSProblem{FT}, trial_points::Vector{Vector{FT}})::IterationOutcome where {FT<:AbstractFloat}

Multi-threaded evaluation of set of trial points. Uses the number of threads that Julia was started with.
"""
function EvaluatePointParallel!(p::DSProblem{FT}, trial_points::Vector{Vector{FT}})::IterationOutcome where {FT<:AbstractFloat}
    #TODO could split into an evaluation function and an update function
    isempty(trial_points) && return Unsuccessful

    #Variables to store the best evaluated points
    feasible_point = isnothing(p.x) ? FT(Inf) * ones(p.N) : p.x
    feasible_cost = isnothing(p.x_cost) ? Threads.Atomic{FT}(FT(Inf)) : Threads.Atomic{FT}(p.x_cost)
    infeasible_point = isnothing(p.i) ? FT(Inf) * ones(p.N) : p.i
    infeasible_cost = isnothing(p.i_cost) ? Threads.Atomic{FT}(FT(Inf)) : Threads.Atomic{FT}(p.i_cost)
    successful_direction = nothing

    #The current minimum hmax value for all collections
    h_min = GetOldHmaxSum(p.constraints)

    updated = Threads.Atomic{Bool}(false)

    #Iterate over all trial points
    Threads.@threads for i=1:length(trial_points)
        point = trial_points[i]

        p.config.opportunistic && updated[] && break

        feasibility = ConstraintEvaluation(p, point)

        # Point violates h_max on at least one constraint collection
        feasibility == StrongInfeasible && continue

        #Point is feasible for relaxed constraints, so evaluate it
        p.status.blackbox_time_total += @elapsed (cost, is_from_cache) = function_evaluation_parallel(p, point)

        #To determine if a point is dominating or improving the combined h_max is needed
        h = GetViolationSumParallel(p, point)

        # Conditions met for a dominant point
        if feasibility == Feasible && cost < feasible_cost[]
            feasible_point = point
            Threads.atomic_xchg!(feasible_cost, cost)
            successful_direction = isnothing(p.status.directions) ? nothing : p.status.directions[i]
            # updated = true
            Threads.atomic_xchg!(updated, true)
        elseif feasibility == WeakInfeasible && h < h_min
            # Conditions met for an improving point (worse cost, but closer to being feasible) or
            # a dominant point (better cost and closer to feasibility)
            # Only record if it offers an improved constraint violation
            infeasible_point = point
            Threads.atomic_xchg!(infeasible_cost, cost)
            successful_direction = isnothing(p.status.directions) ? nothing : p.status.directions[i]
            h_min = h
            Threads.atomic_xchg!(updated, true)
        end

        p.full_output && OutputPointEvaluation(i, point, cost, h, is_from_cache, Threads.threadid())
    end

    result = Unsuccessful

    incum_i_cost = isnothing(p.i_cost) ? FT(Inf) : p.i_cost
    incum_x_cost = isnothing(p.x_cost) ? FT(Inf) : p.x_cost

    p.status.directions = nothing

    # Dominates if there is a feasible improvement, or an infeasible point with
    # reduced violation (determined previously) as well as a cost lower than any
    # yet tested (feasible and infeasible)
    if feasible_cost[] < incum_x_cost
        result = Dominating
    elseif infeasible_cost[] < incum_i_cost && h_min < GetOldHmaxSum(p.constraints)
        result = Dominating
    elseif h_min < GetOldHmaxSum(p.constraints)
        result = Improving
    end

    if feasible_cost[] < incum_x_cost
        p.x = feasible_point
        p.x_cost = feasible_cost[]
    end

    if infeasible_cost[] < incum_i_cost || h_min < GetOldHmaxSum(p.constraints)
        p.i = infeasible_point
        p.i_cost = infeasible_cost[]
    end

    p.status.success_direction = successful_direction

    UpdateConstraints(p.constraints, h_min, result, p.x, p.i)

    return result
end

#Wrapper for matching to empty trial point arrays
(EvaluatePoint!(p::DSProblem{T}, trial_points::Vector)::IterationOutcome) where T =
    EvaluatePoint!(p, convert(Vector{Vector{T}}, trial_points))

"""
	function_evaluation(p::DSProblem{T}, trial_point::Vector{T})::(T, Bool) where T

Evaluate a single trial point with the objective function of `p`.

By default calls the function with the trial point and returns the result. Override to
provide custom evaluation behaviour.
"""
function function_evaluation(p::DSProblem{T},
                             trial_point::Vector{T})::Tuple{T, Bool} where T
    if CacheQuery(p, trial_point)
        p.status.cache_hits += 1
        return (CacheGet(p, trial_point), true)
    end
    cost = round(p.objective(trial_point), digits=p.config.cost_digits)
    p.status.function_evaluations += 1
    CachePush(p, trial_point, cost)
	return (cost, false)
end

"""
    function_evaluation_parallel(p::DSProblem{T}, trial_point::Vector{T})::Tuple{T, Bool} where T

Evaluate a single trial point with the objective function of `p` using multiple threads.
"""
function function_evaluation_parallel(p::DSProblem{T}, trial_point::Vector{T})::Tuple{T, Bool} where T
    if CacheQueryParallel(p, trial_point)
        p.status.cache_hits += 1
        return (CacheGetParallel(p, trial_point), true)
    end
    cost = p.objective(trial_point)
    p.status.function_evaluations += 1
    CachePushParallel(p, trial_point, cost)
	return (cost, false)
end

function _check_initial_point(p::DSProblem{T}) where T
    for i=1:p.N
        if p.granularity[i] > 0 && (p.user_initial_point[i] / p.granularity[i]) % 1 != 0
            error("Initial value of variable with index $i is not an integer multiple of its granularity.")
        end
    end
end


text:
module TConv2DM
    include("../tools/glorotuniform.jl")
    using .GlorotUniform
    using LoopVectorization

    mutable struct TConv2D
        save_layer::Any
        load_layer::Any
        activate::Any
        initialize::Any
        update::Any

        input_shape::Tuple
        output_shape::Tuple
        kernel_size::Tuple{Int64, Int64}
        activation_function::Module

        filters::Array{Float32}
        biases::Array{Float32}

        padding::Int64
        filter::Int64
        strides::Tuple{Int64, Int64}

        ∇weights::Array{Float32}
        Vdw::Array{Float32}
        Sdw::Array{Float32}
        ∇biases_c::Array{Float32}
        Vdb::Array{Float32}
        Sdb::Array{Float32}

        padding_input::Array{Float32}
        value::Array{Float32}
        output::Array{Float32}

        ∇biases::Array{Float32}
        pre_δ::Array{Float32}
        δ::Array{Float32}

        function TConv2D(;input_shape::Tuple, filter::Int64, padding::Int64=0, kernel_size::Tuple{Int64, Int64}, strides::Tuple{Int64, Int64}=(1,1), activation_function::Module, randomization::Bool=true)
            conv_num_per_col = (input_shape[1]-1)*strides[1]+2*padding-kernel_size[1]+2
            conv_num_per_row = (input_shape[2]-1)*strides[2]+2*padding-kernel_size[2]+2
            output_shape = (conv_num_per_col, conv_num_per_row, filter)

            input_size = 1
            for i in input_shape
                input_size *= i
            end

            filter_size = (filter, input_shape[3], kernel_size[1], kernel_size[2])
            filters = !randomization ? zeros(Float32, filter_size) : Array{Float32}(GUMatrix(input_size, conv_num_per_row*conv_num_per_col*filter, filter_size))
            biases = !randomization ? zeros(Float32, filter) : Array{Float32}(GUMatrix(input_size, conv_num_per_row*conv_num_per_col*filter, filter))

            ∇weights = zeros(Float32, filter_size)
            Vdw = zeros(Float32, filter_size)
            Sdw = zeros(Float32, filter_size)
            ∇biases_c = zeros(Float32, filter)
            Vdb = zeros(Float32, filter)
            Sdb = zeros(Float32, filter)
            new(save_TConv2D, load_TConv2D, activate_TConv2D, init_TConv2D, update_TConv2D, input_shape, output_shape, kernel_size, activation_function, filters, biases, padding, filter, strides, ∇weights, Vdw, Sdw, ∇biases_c, Vdb, Sdb)
        end
    end

    function init_TConv2D(layer::TConv2D, mini_batch::Int64)
        layer.padding_input = zeros(Float32, ((layer.input_shape[1]-1)*layer.strides[1]+2*layer.padding+1, (layer.input_shape[2]-1)*layer.strides[2]+2*layer.padding+1, layer.input_shape[3], mini_batch))
        layer.value = zeros(Float32, layer.output_shape..., mini_batch)
        layer.output = zeros(Float32, layer.output_shape..., mini_batch)

        layer.∇biases = zeros(Float32, layer.output_shape..., mini_batch)
        layer.pre_δ = zeros(Float32, ((layer.input_shape[1]-1)*layer.strides[1]+2*layer.padding+1, (layer.input_shape[2]-1)*layer.strides[2]+2*layer.padding+1, layer.input_shape[3], mini_batch))
        layer.δ = zeros(Float32, layer.input_shape..., mini_batch)
    end

    function activate_TConv2D(layer::TConv2D, input::Array{Float32})
        x, y = layer.strides
        @avxt for i in axes(input, 1), j in axes(input, 2), c in axes(input, 3), b in axes(input, 4)
            layer.padding_input[layer.padding+(i-1)*x+1, layer.padding+(j-1)*y+1, c, b] = input[i, j, c, b]
        end

        @avxt for i in axes(layer.value, 1), j in axes(layer.value, 2), f in axes(layer.value, 3), b in axes(layer.value, 4)
            s = 0.0f0
            for k₁ in axes(layer.filters, 3), k₂ in axes(layer.filters, 4), c in axes(layer.filters, 2)
                s += layer.padding_input[i-1+k₁, j-1+k₂, c, b] * layer.filters[f, c, k₁, k₂]
            end
            layer.value[i, j, f, b] = s + layer.biases[f]
        end

        layer.activation_function.func!(layer.output, layer.value)
    end

    function update_TConv2D(layer::TConv2D, optimizer::String, input::Array{Float32}, δₗ₊₁::Array{Float32}, α::Float64, parameters::Tuple, direction::Int64=1)
        layer.activation_function.get_∇biases!(layer.∇biases, layer.value, δₗ₊₁)
        @avxt layer.pre_δ .= 0.0f0

        @avxt for i in axes(layer.∇biases, 1), j in axes(layer.∇biases, 2), f in axes(layer.∇biases, 3), b in axes(layer.∇biases, 4)
            for k₁ in axes(layer.filters, 3), k₂ in axes(layer.filters, 4), c in axes(layer.filters, 2)
                layer.pre_δ[i-1+k₁, j-1+k₂, c, b] += layer.filters[f, c, k₁, k₂]*layer.∇biases[i, j, f, b]
            end
        end

        x, y = layer.strides
        @avxt for i in axes(layer.δ, 1), j in axes(layer.δ, 2), c in axes(layer.δ, 3), b in axes(layer.δ, 4)
            layer.δ[i, j, c, b] = layer.pre_δ[layer.padding+(i-1)*x+1, layer.padding+(j-1)*y+1, c, b]
        end

        if optimizer=="SGD"
            @avxt for i in axes(layer.∇biases, 1), j in axes(layer.∇biases, 2), f in axes(layer.∇biases, 3)
                for k₁ in axes(layer.filters, 3), k₂ in axes(layer.filters, 4), c in axes(layer.filters, 2)
                    layer.filters[f, c, k₁, k₂] -= α*layer.padding_input[(i-1)*x+k₁, (j-1)*y+k₂, c, 1]*layer.∇biases[i, j, f, 1]*direction
                end
                layer.biases[f] -= α*layer.∇biases[i, j, f, 1]*direction
            end

        elseif optimizer=="Minibatch_GD"
            batch_size = size(layer.∇biases, 4)
            @avxt for i in axes(layer.∇biases, 1), j in axes(layer.∇biases, 2), f in axes(layer.∇biases, 3), b in axes(layer.∇biases, 4)
                for k₁ in axes(layer.filters, 3), k₂ in axes(layer.filters, 4), c in axes(layer.filters, 2)
                    layer.filters[f, c, k₁, k₂] -= α*layer.padding_input[(i-1)*x+k₁, (j-1)*y+k₂, c, b]*layer.∇biases[i, j, f, b]/batch_size*direction
                end
                layer.biases[f] -= α*layer.∇biases[i, j, f, b]/batch_size*direction
            end

        elseif optimizer=="Adam"
            t = parameters[1]
            β₁ = parameters[2]
            β₂ = parameters[3]
            ϵ = parameters[4]

            @avxt layer.∇weights .= 0.0f0
            @avxt layer.∇biases_c .= 0.0f0

            @avxt for i in axes(layer.∇biases, 1), j in axes(layer.∇biases, 2), f in axes(layer.∇biases, 3)
                for k₁ in axes(layer.filters, 3), k₂ in axes(layer.filters, 4), c in axes(layer.filters, 2)
                    layer.∇weights[f, c, k₁, k₂] += layer.padding_input[(i-1)*x+k₁, (j-1)*y+k₂, c, 1]*layer.∇biases[i, j, f, 1]
                end
                layer.∇biases_c[f] += layer.∇biases[i, j, f, 1]
            end

            @avxt for i in eachindex(layer.filters)
                layer.Vdw[i] = β₁*layer.Vdw[i] + (1-β₁)*layer.∇weights[i]
                layer.Sdw[i] = β₂*layer.Sdw[i] + (1-β₂)*layer.∇weights[i]^2
                layer.filters[i] -= α*(layer.Vdw[i]/(1-β₁^t))/(sqrt(layer.Sdw[i]/(1-β₂^t))+ϵ)*direction
            end
            @avxt for i in eachindex(layer.biases)
                layer.Vdb[i] = β₁*layer.Vdb[i] + (1-β₁)*layer.∇biases_c[i]
                layer.Sdb[i] = β₂*layer.Sdb[i] + (1-β₂)*layer.∇biases_c[i]^2
                layer.biases[i] -= α*(layer.Vdb[i]/(1-β₁^t))/(sqrt(layer.Sdb[i]/(1-β₂^t))+ϵ)*direction
            end

        elseif optimizer=="AdaBelief"
            t = parameters[1]
            β₁ = parameters[2]
            β₂ = parameters[3]
            ϵ = parameters[4]

            @avxt layer.∇weights .= 0.0f0
            @avxt layer.∇biases_c .= 0.0f0

            @avxt for i in axes(layer.∇biases, 1), j in axes(layer.∇biases, 2), f in axes(layer.∇biases, 3)
                for k₁ in axes(layer.filters, 3), k₂ in axes(layer.filters, 4), c in axes(layer.filters, 2)
                    layer.∇weights[f, c, k₁, k₂] += layer.padding_input[(i-1)*x+k₁, (j-1)*y+k₂, c, 1]*layer.∇biases[i, j, f, 1]
                end
                layer.∇biases_c[f] += layer.∇biases[i, j, f, 1]
            end

            @avxt for i in eachindex(layer.filters)
                layer.Vdw[i] = β₁*layer.Vdw[i] + (1-β₁)*layer.∇weights[i]
                layer.Sdw[i] = β₂*layer.Sdw[i] + (1-β₂)*(layer.∇weights[i]-layer.Vdw[i])^2
                layer.filters[i] -= α*(layer.Vdw[i]/(1-β₁^t))/(sqrt(layer.Sdw[i]/(1-β₂^t))+ϵ)*direction
            end
            @avxt for i in eachindex(layer.biases)
                layer.Vdb[i] = β₁*layer.Vdb[i] + (1-β₁)*layer.∇biases_c[i]
                layer.Sdb[i] = β₂*layer.Sdb[i] + (1-β₂)*(layer.∇biases_c[i]-layer.Vdb[i])^2
                layer.biases[i] -= α*(layer.Vdb[i]/(1-β₁^t))/(sqrt(layer.Sdb[i]/(1-β₂^t))+ϵ)*direction
            end
        end
    end

    function save_TConv2D(layer::TConv2D, file::Any, id::String)
        write(file, id, "TConv2D")
        write(file, id*"input_shape", collect(layer.input_shape))
        write(file, id*"filter", layer.filter)
        write(file, id*"padding", layer.padding)
        write(file, id*"kernel_size", collect(layer.kernel_size))
        write(file, id*"strides", collect(layer.strides))
        write(file, id*"activation_function", layer.activation_function.get_name())

        write(file, id*"filters", layer.filters)
        write(file, id*"biases", layer.biases)
    end

    function load_TConv2D(layer::TConv2D, file::Any, id::String)
        layer.filters = read(file, id*"filters")
        layer.biases = read(file, id*"biases")
    end

    function get_args(file::Any, id::String)
        input_shape = tuple(read(file, id*"input_shape")...)
        filter = read(file, id*"filter")
        padding = read(file, id*"padding")
        kernel_size = tuple(read(file, id*"kernel_size")...)
        strides = tuple(read(file, id*"strides")...)
        return (input_shape=input_shape, filter=filter, padding=padding, kernel_size=kernel_size, strides=strides)
    end
end

ext: jl


text:
using EqualityConstrainedMinimization
using Test
using LinearAlgebra

#Test objective for constrained problem (see Boyd and Vandenberghe problem 10.15)
function Test_Objective_Con(_x)
    _n = length(_x)
    sum = 0

    for i = 1:_n
        sum += _x[i] * log(_x[i])
    end

    return sum
end

#Generates ECM problem with feasible starting point
function Generate_ECM(_n, _p, A_bound)
    _x = rand(_n)
    _A = rand(_p, _n) * A_bound

    it = 0
    #Ensure A has full rank
    while rank(_A) != max(_n, _p) && it < 100 #Cap generation at 100 tries. The odds of needing this many are very low
        _A = rand(_p, _n) * A_bound
        it += 1
    end

    _b = _A * _x

    return _A, _x, _b
end

l = 3
w = 2
mtx = [0.52 0.96 0.59; 0.42 0.37 0.25]
feasX = [0.2; 0.3; 0.4]
x = ones(l)
vect = mtx*feasX

n = 100
p = 30
ECM = Generate_ECM(n, p, 1)

f(_x) = Test_Objective_Con(_x)

@testset "EqualityConstrainedMinimization.jl" begin
    @test Solve_ECM(f, x, verbose = false) ≈ -1.1036383232 rtol = 1e-8
    @test Solve_ECM(f, x, A = mtx, b = vect, verbose = false) ≈ -1.0587361852 rtol = 1e-8
    @test Solve_ECM(f, feasX, A = mtx, verbose = false) ≈ -1.0587361852 rtol = 1e-8
    @test Solve_ECM(f, feasX, A = mtx, b = vect, verbose = false) ≈ -1.0587361852 rtol = 1e-8

    @test Solve_ECM(f, ECM[2], A = ECM[1], verbose = false) ≈ Solve_ECM(f, ones(n), A = ECM[1], b = ECM[3], verbose = false) rtol = 1e-8
end

max_stars_repo_name: ihawn/EqualityConstrainedMinimization

max_issues_repo_name: ihawn/EqualityConstrainedMinimization

max_forks_repo_name: ihawn/EqualityConstrainedMinimization


text:
#-----------------------------------------------------------------#
# This routine computes the RHS forcing of a simple two-rate ODE
# Written by F.X. Giraldo / P.R. Mugg on 9/19/2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey, CA 93943-5216
#------------------------------------------------------------------#
function rhs_function(q0,c,time)

    # RHS components
    f = c*q0;        #fast component
    g = sin(time);   #slow component

    # Construct RHS
    rhs = f + g;

    return (rhs,f,g)
end #end of rhs_function


text:
using UnicodePlots
using Test
using SparseArrays: sprand
import Dates: Date, Day

# write your own tests here
include("tst_plots.jl")


text:
using PyPlot

function draw3Graph()

	file = open("../Experimentation/"*"ConvexeBase", "r")

	moyTimeConvexeBase = zeros(Float64, 11)
	moyNbIterConvexeBase = zeros(Float64, 11)
	nbInstConvexeBase = zeros(Int, 11)

	x = readline(file)
	while x != ""
		informations = split(x, ";")
		nbVar = parse(Int, informations[1])
		numInst = parse(Int, informations[2])
		timeTaken = parse(Float64, informations[3])
		nbIter = parse(Int, informations[4])

		moyTimeConvexeBase[Int(nbVar/10 - 4)] += timeTaken
		moyNbIterConvexeBase[Int(nbVar/10 - 4)] += nbIter
		nbInstConvexeBase[Int(nbVar/10 - 4)] += 1

		x = readline(file)
	end
	moyTimeConvexeBase ./= nbInstConvexeBase
	moyNbIterConvexeBase ./= nbInstConvexeBase

	close(file)

	file = open("../Experimentation/LinearBase", "r")

	moyTimeLinearBase = zeros(Float64, 11)
	moyNbIterLinearBase = zeros(Float64, 11)
	nbInstLinearBase = zeros(Int, 11)

	x = readline(file)
	while x != ""
		informations = split(x, ";")
		nbVar = parse(Int, informations[1])
		numInst = parse(Int, informations[2])
		timeTaken = parse(Float64, informations[3])
		nbIter = parse(Int, informations[4])

		moyTimeLinearBase[Int(nbVar/10 - 4)] += timeTaken
		moyNbIterLinearBase[Int(nbVar/10 - 4)] += nbIter
		nbInstLinearBase[Int(nbVar/10 - 4)] += 1

		x = readline(file)
	end
	moyTimeLinearBase ./= nbInstLinearBase
	moyNbIterLinearBase ./= nbInstLinearBase

	close(file)

	file = open("../Experimentation/Combo", "r")

	moyTimeCombo = zeros(Float64, 11)
	moyNbIterCombo = zeros(Float64, 11)
	nbInstCombo = zeros(Int, 11)

	x = readline(file)
	while x != ""
		informations = split(x, ";")
		nbVar = parse(Int, informations[1])
		numInst = parse(Int, informations[2])
		timeTaken = parse(Float64, informations[3])
		nbIter = parse(Int, informations[4])

		moyTimeCombo[Int(nbVar/10 - 4)] += timeTaken
		moyNbIterCombo[Int(nbVar/10 - 4)] += nbIter
		nbInstCombo[Int(nbVar/10 - 4)] += 1

		x = readline(file)
	end
	moyTimeCombo ./= nbInstCombo
	moyNbIterCombo ./= nbInstCombo

	close(file)

	file = open("../Experimentation/ComboHeuristic", "r")

	moyTimeComboHeuristic = zeros(Float64, 11)
	moyNbIterComboHeuristic = zeros(Float64, 11)
	nbInstComboHeuristic = zeros(Int, 11)

	x = readline(file)
	while x != ""
		informations = split(x, ";")
		nbVar = parse(Int, informations[1])
		numInst = parse(Int, informations[2])
		timeTaken = parse(Float64, informations[3])
		nbIter = parse(Int, informations[4])

		moyTimeComboHeuristic[Int(nbVar/10 - 4)] += timeTaken
		moyNbIterComboHeuristic[Int(nbVar/10 - 4)] += nbIter
		nbInstComboHeuristic[Int(nbVar/10 - 4)] += 1

		x = readline(file)
	end
	moyTimeComboHeuristic ./= nbInstComboHeuristic
	moyNbIterComboHeuristic ./= nbInstComboHeuristic

	close(file)

	file = open("../Experimentation/LinearHeuristic", "r")

	moyTimeLinearHeuristic = zeros(Float64, 11)
	moyNbIterLinearHeuristic = zeros(Float64, 11)
	nbInstLinearHeuristic = zeros(Int, 11)

	x = readline(file)
	while x != ""
		informations = split(x, ";")
		nbVar = parse(Int, informations[1])
		numInst = parse(Int, informations[2])
		timeTaken = parse(Float64, informations[3])
		nbIter = parse(Int, informations[4])

		moyTimeLinearHeuristic[Int(nbVar/10 - 4)] += timeTaken
		moyNbIterLinearHeuristic[Int(nbVar/10 - 4)] += nbIter
		nbInstLinearHeuristic[Int(nbVar/10 - 4)] += 1

		x = readline(file)
	end
	moyTimeLinearHeuristic ./= nbInstLinearHeuristic
	moyNbIterLinearHeuristic ./= nbInstLinearHeuristic

	close(file)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeConvexeBase[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeConvexeBase[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Temps de résolution moyen (en secondes)")
	ax.grid(true)
	fig.savefig("SaveFig/moyTimeConvexeBase.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeConvexeBase[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyNbIterConvexeBase[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Nombre de sous-problème moyen étudié")
	ax.grid(true)
	fig.savefig("SaveFig/moyNbIterConvexeBase.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeConvexeBase[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeConvexeBase[Int.((solX ./ 10) .- 4)] ./ moyNbIterConvexeBase[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Temps de résolution moyen d'un sous-problème (en secondes)")
	ax.grid(true)
	fig.savefig("SaveFig/moyTimeByIterConvexeBase.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeLinearBase[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeLinearBase[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Temps de résolution moyen (en secondes)")
	ax.grid(true)
	fig.savefig("SaveFig/moyTimeLinearBase.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeLinearBase[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyNbIterLinearBase[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Nombre de sous-problème moyen étudié")
	ax.grid(true)
	fig.savefig("SaveFig/moyNbIterLinearBase.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeLinearBase[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeLinearBase[Int.((solX ./ 10) .- 4)] ./ moyNbIterLinearBase[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Temps de résolution moyen d'un sous-problème (en secondes)")
	ax.grid(true)
	fig.savefig("SaveFig/moyTimeByIterLinearBase.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeCombo[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeCombo[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Temps de résolution moyen (en secondes)")
	ax.grid(true)
	fig.savefig("SaveFig/moyTimeCombo.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeCombo[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyNbIterCombo[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Nombre de sous-problème moyen étudié")
	ax.grid(true)
	fig.savefig("SaveFig/moyNbIterCombo.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeCombo[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeCombo[Int.((solX ./ 10) .- 4)] ./ moyNbIterCombo[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Temps de résolution moyen d'un sous-problème (en secondes)")
	ax.grid(true)
	fig.savefig("SaveFig/moyTimeByIterCombo.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeComboHeuristic[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeComboHeuristic[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Temps de résolution moyen (en secondes)")
	ax.grid(true)
	fig.savefig("SaveFig/moyTimeComboHeuristic.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeComboHeuristic[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyNbIterComboHeuristic[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Nombre de sous-problème moyen étudié")
	ax.grid(true)
	fig.savefig("SaveFig/moyNbIterComboHeuristic.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeComboHeuristic[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeComboHeuristic[Int.((solX ./ 10) .- 4)] ./ moyNbIterComboHeuristic[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Temps de résolution moyen d'un sous-problème (en secondes)")
	ax.grid(true)
	fig.savefig("SaveFig/moyTimeByIterComboHeuristic.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeLinearHeuristic[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeLinearHeuristic[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Temps de résolution moyen (en secondes)")
	ax.grid(true)
	fig.savefig("SaveFig/moyTimeLinearHeuristic.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeLinearHeuristic[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyNbIterLinearHeuristic[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Nombre de sous-problème moyen étudié")
	ax.grid(true)
	fig.savefig("SaveFig/moyNbIterLinearHeuristic.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeLinearHeuristic[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeLinearHeuristic[Int.((solX ./ 10) .- 4)] ./ moyNbIterLinearHeuristic[Int.((solX ./ 10) .- 4)]

	ax.plot(solX, solY, "b.-")
	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Temps de résolution moyen d'un sous-problème (en secondes)")
	ax.grid(true)
	fig.savefig("SaveFig/moyTimeByIterLinearHeuristic.png")
	close(fig)

	###############################################################################""

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeConvexeBase[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeConvexeBase[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Convexe JuMP")

	solX = filter(x->moyTimeLinearBase[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeLinearBase[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Relaxation Linéaire")

	solX = filter(x->moyTimeCombo[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeCombo[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Combo")

	solX = filter(x->moyTimeComboHeuristic[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeComboHeuristic[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Combo + Heuristique Primale")

	solX = filter(x->moyTimeLinearHeuristic[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeLinearHeuristic[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Relaxation linéaire + Heuristique Primale")

	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Temps de résolution moyen (en secondes)")
	ax.grid(true)
	ax.legend()
	fig.savefig("SaveFig/moyTimeGeneral.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyNbIterConvexeBase[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyNbIterConvexeBase[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Convexe JuMP")

	solX = filter(x->moyNbIterLinearBase[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyNbIterLinearBase[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Relaxation Linéaire")

	solX = filter(x->moyNbIterCombo[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyNbIterCombo[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Combo")

	solX = filter(x->moyNbIterComboHeuristic[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyNbIterComboHeuristic[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Combo + Heuristique Primale")

	solX = filter(x->moyNbIterLinearHeuristic[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyNbIterLinearHeuristic[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Relaxation linéaire + Heuristique Primale")

	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Nombre de sous-problème étudié")
	ax.grid(true)
	ax.legend()
	fig.savefig("SaveFig/moyNbIterGeneral.png")
	close(fig)

	fig = figure()
	ax = fig.add_subplot(111)

	solX = filter(x->moyTimeConvexeBase[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeConvexeBase[Int.((solX ./ 10) .- 4)] ./ moyNbIterConvexeBase[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Convexe JuMP")

	solX = filter(x->moyTimeLinearBase[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeLinearBase[Int.((solX ./ 10) .- 4)] ./ moyNbIterLinearBase[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Relaxation Linéaire")

	solX = filter(x->moyTimeCombo[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeCombo[Int.((solX ./ 10) .- 4)] ./ moyNbIterCombo[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Combo")

	solX = filter(x->moyTimeComboHeuristic[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeComboHeuristic[Int.((solX ./ 10) .- 4)] ./ moyNbIterComboHeuristic[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Combo + Heuristique Primale")

	solX = filter(x->moyTimeLinearHeuristic[Int(x/10 - 4)]!=0, 50:10:150)
	solY = moyTimeLinearHeuristic[Int.((solX ./ 10) .- 4)] ./ moyNbIterLinearHeuristic[Int.((solX ./ 10) .- 4)]
	ax.plot(solX, solY, ".-", label = "Relaxation linéaire + Heuristique Primale")

	ax.set_xlabel("Nombre d'objets")
	ax.set_ylabel("Temps de résolution moyen d'un sous-problème (en secondes)")
	ax.grid(true)
	ax.legend()
	fig.savefig("SaveFig/moyTimeByIterGeneral.png")
	close(fig)

end

max_stars_repo_name: LucasBaussay/BranchAndBoundBiObjKp

max_issues_repo_name: LucasBaussay/BranchAndBoundBiObjKp

max_forks_repo_name: LucasBaussay/BranchAndBoundBiObjKp


text:
using JuMP, Gurobi, Distributions, StatsFuns

const gurobi_env = Gurobi.Env()

# DIMENSIONS
const numVariables = 100
const numJCC = 100
const conRHS = numVariables^2

# CASES
# number of samples, replicates for scenario approximations
const NumSamples = ceil.(Int64,logspace(1,log10(50000.0),50))
const NumReplicates = 20

# FILENAMES
const objFile = "objectiveValue.txt"
const solnFile = "solution.txt"
const riskFile = "riskLevel.txt"
const riskTimeFile = "riskTime.txt"
const solnTimeFile = "solutionTime.txt"
const checkingTimeFile = "checkingTime.txt"
const sortingTimeFile = "sortingTime.txt"
const callbackFile = "numCallbacks.txt"
const numConstraintsFile = "numConstraintsEnforced.txt"

xi_mean = zeros(Float64,numVariables)
xi_covariance = 0.5*ones(Float64,numJCC,numJCC)

for j = 1:numVariables
	xi_mean[j] = j*1.0/numVariables
end

for i = 1:numJCC
	xi_covariance[i,i] = 1.0
end

xi_cov_chol_tmp = cholfact(xi_covariance)
xi_cov_chol = xi_cov_chol_tmp[:L]

xi_mean_mat = zeros(Float64,numJCC,numVariables)
for j = 1:numVariables
	xi_mean_mat[:,j] = xi_mean[j]*ones(numJCC)
end

# tailored implementation of multivariate normal distribution
function myMvNormal(numSamples::Int64)

	xi = zeros(Float64,numJCC,numVariables,numSamples)

	xi_tmp = rand(Normal(0.0,1.0),numJCC,numVariables,numSamples)
	for samp = 1:numSamples
		xi[:,:,samp] = xi_cov_chol*xi_tmp[:,:,samp] + xi_mean_mat
	end

	return xi
end

# generate random samples from multivariate normal distribution
function generateRandomSamples(numSamples::Int64)

	xi = myMvNormal(numSamples)

	return xi
end

srand(1234)

# generate samples to estimate actual risk level of solution
const numAnalyticalSamples = 20000
const xi_analytical = generateRandomSamples(numAnalyticalSamples)

srand()

# compute negative probability of all constraints being satisfied
# for given scenarios of demands d and a given decision vector x
function computeRiskLevel(x::Array{Float64},reliabilityLevel::Float64=1E-06)

	const numSamples = size(xi_analytical,3)
	simpleEst::Bool = false

	numViolated::Int64 = checkScenarioConstraints(x,xi_analytical)

	riskLevel::Float64 = Inf

	if(simpleEst)
		riskLevel = numViolated*1.0/numSamples
		return riskLevel
	else
		if(numViolated == 0)
			riskLevel = -log(reliabilityLevel)/numSamples
			return riskLevel
		end

		 # termination criterion for bisection
		gamma_tolerance::Float64 = 0.1/numSamples

		# use tail bounds on binomial distribution to estimate lower and upper bounds
		gamma_low::Float64 = (numViolated - sqrt(-numSamples*log(reliabilityLevel)/2.0))/numSamples
		if(gamma_low < 0.0)
			gamma_low = 0.0
		end
		gamma_up::Float64 = (numViolated + sqrt(-numSamples*log(reliabilityLevel)/2.0))/numSamples
		if(gamma_up > 1.0)
			gamma_up = 1.0
		end

		gamma::Float64 = 0.0
		func_value::Float64 = 0.0

		# use bisection to solve the nonlinear equation
		while (gamma_up - gamma_low > gamma_tolerance)

			gamma = (gamma_low + gamma_up)/2.0
			func_value = -reliabilityLevel
			for i = 1:numViolated
				tmp_func = i*log(gamma) + (numSamples-i)*log(1-gamma) + lfact(numSamples) - lfact(numSamples-i) - lfact(i)
				func_value += exp(tmp_func)
				if(func_value > 0)
					break
				end
			end
			if(func_value > 0)
				gamma_low = (gamma_low + gamma_up)/2.0
			else
				gamma_up = (gamma_low + gamma_up)/2.0
			end

		end

		riskLevel = gamma_up
		return riskLevel
	end

end

# determines if all constraints are satisfied for a given single scenario
# returns one if scenario constraints are all satisfied, zero otherwise
function checkScenarioConstraints(x::Array{Float64},xi::Array{Float64})

	numSamples::Int64 = size(xi,3)

	numViolated::Int64 = 0
	for iter = 1:numSamples
		for i = 1:numJCC
			con::Float64 = evaluateConstraint(i,x,xi[i,:,iter])
			if(con > 0)
				numViolated += 1
				break
			end
		end
	end

	return numViolated
end

# evaluate the constraint values for given decision vector x,
# realization of the random variables xi, and constraint scalings
function evaluateConstraint(index::Int64,x::Array{Float64},xi::Array{Float64})

	con::Float64 = norm(xi.*x)^2 - conRHS

	return con
end

# determines if all constraints are satisfied for a given single scenario
# returns one if scenario constraints are all satisfied, zero otherwise
function getScenarioConstraintViolations(x::Array{Float64},xi::SubArray{Float64,2,Array{Float64,3},Tuple{Int64,Base.Slice{Base.OneTo{Int64}},Base.Slice{Base.OneTo{Int64}}},true},pickFirstSetOfConstraints::Bool,maxNumConstraintsPerIteration::Int64)

	const numSamples = size(xi,2)
	numViolated::Int64 = 0
	constraintViolations = Float64[]
	constraintViolationIndices = Int64[]
	constraintViolationTolerance::Float64 = 1E-06

	for iter = 1:numSamples
		@views con::Float64 = norm(xi[:,iter].*x) - sqrt(conRHS)

		if(con >= constraintViolationTolerance)
			numViolated += 1
			push!(constraintViolations,con)
			push!(constraintViolationIndices,iter)
			if(pickFirstSetOfConstraints && numViolated == maxNumConstraintsPerIteration)
				break
			end
		end
	end

	return numViolated, constraintViolations, constraintViolationIndices
end

# write solution to file
function initializeFiles(details_file::String)

	writeMode::String = "w"

	open(details_file, writeMode) do f
		write(f,"Model: $modelName \n")
		write(f,"numVariables: $numVariables \n")
		write(f,"numJCC: $numJCC \n")
		write(f,"conRHS: $conRHS \n")
		write(f,"NumSamples: $NumSamples \n")
		write(f,"NumReplicates: $NumReplicates \n")
		write(f,"Method: Scenario approximations \n")
	end

end

# write solution to file
function writeSolutionToFile(dirName::String, status::Symbol, objvar_opt::Float64, x_opt::Array{Float64}, riskLevel::Float64, solutionTime::Float64, checkingTime::Float64, sortingTime::Float64, numCallbacks::Int64, numConstraintsAdded::Int64, riskTime::Float64)

	if(status == Symbol("Optimal"))
		writeMode::String = "a"

		obj_file::String = dirName * objFile
		open(obj_file, writeMode) do f
			write(f,"$objvar_opt \n")
		end

		solution_file::String = dirName * solnFile
		open(solution_file, writeMode) do f
			for i = 1:numVariables
				write(f,"$(x_opt[i]) ")
			end
			write(f,"\n")
		end

		risk_file::String = dirName * riskFile
		open(risk_file, writeMode) do f
			write(f,"$riskLevel \n")
		end

		risktime_file::String = dirName * riskTimeFile
		open(risktime_file, writeMode) do f
			write(f,"$riskTime \n")
		end

		solntime_file::String = dirName * solnTimeFile
		open(solntime_file, writeMode) do f
			write(f,"$solutionTime \n")
		end

		checktime_file::String = dirName * checkingTimeFile
		open(checktime_file, writeMode) do f
			write(f,"$checkingTime \n")
		end

		sorttime_file::String = dirName * sortingTimeFile
		open(sorttime_file, writeMode) do f
			write(f,"$sortingTime \n")
		end

		callback_file::String = dirName * callbackFile
		open(callback_file, writeMode) do f
			write(f,"$numCallbacks \n")
		end

		constraints_file::String = dirName * numConstraintsFile
		open(constraints_file, writeMode) do f
			write(f,"$numConstraintsAdded \n")
		end
	end

end

# solve the scenario approximation model
function solveScenarioApproximationModel(numScenarios::Int64,dirName::String)

	tic()

	# sample the random variables
	xi = generateRandomSamples(numScenarios)

	maxNumConstraintsPerIteration::Int64 = 10
	numConstraintsAdded::Int64 = 0

	pickFirstSetOfConstraints::Bool = false

	sortingTime::Float64 = 0.0
	checkingTime::Float64 = 0.0
	numCallbacks::Int64 = 0

	objvar_opt::Float64 = Inf
	y_opt = zeros(Float64,numVariables)
	status::Symbol = Symbol("Error")

	numViolations::Int64 = numScenarios

	constraintIndices = [Int64[] for i=1:numJCC]

	while(numViolations > 0)

		numCallbacks += 1
		numViolations = 0

		numConstraintsConsidered = zeros(Int64,numJCC)
		for i = 1:numJCC
			numConstraintsConsidered[i] = size(constraintIndices[i])[1]
		end

		mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0))

		@variable(mod, objvar)
		@variable(mod, 0 <= y[1:numVariables] <= 10^6)

		@objective(mod, Min, objvar)
		@constraint(mod, -sum(y[j] for j = 1:numVariables) <= objvar)
		@constraint(mod, [i=1:numJCC, scen=1:numConstraintsConsidered[i]], norm(xi[i,:,constraintIndices[i][scen]].*y) <= sqrt(conRHS))

		status = solve(mod)
		objvar_opt = getvalue(objvar)
		y_opt = getvalue(y)

		for i = 1:numJCC

			numViolatedConstraints::Int64 = 0
			constraintViolations = Float64[]
			constraintViolationIndices = Int64[]

			tic()

			@views numViolatedConstraints, constraintViolations, constraintViolationIndices = getScenarioConstraintViolations(y_opt,xi[i,:,:],pickFirstSetOfConstraints,maxNumConstraintsPerIteration)

			checkingTime += toq()

			numConstraintsEnforced::Int64 = min(size(constraintViolations,1),maxNumConstraintsPerIteration)

			tic()
			relevantIndices = 1:1:numConstraintsEnforced
			if(!pickFirstSetOfConstraints)
				relevantIndices = sortperm(constraintViolations, rev=true)
			end
			sortingTime += toq()

			for iter = 1:numConstraintsEnforced
				push!(constraintIndices[i],constraintViolationIndices[relevantIndices[iter]])
			end

			numConstraintsAdded += numConstraintsEnforced

			numViolations += numViolatedConstraints
		end

	end

	tic()

	riskLevel::Float64 = computeRiskLevel(y_opt)

	riskTime = toq()

	solutionTime::Float64 = toq()

	@printf "  Objective: %2.5f,  Risklevel: %.6f,  Time: %4.2f,  NumCallbacks: %d,  NumConstraints: %d \n" objvar_opt riskLevel solutionTime numCallbacks numConstraintsAdded

	writeSolutionToFile(dirName, status, objvar_opt, y_opt, riskLevel, solutionTime, checkingTime, sortingTime, numCallbacks, numConstraintsAdded, riskTime)

	return status, objvar_opt, y_opt, riskLevel, solutionTime, checkingTime, sortingTime, numCallbacks, numConstraintsAdded, riskTime
end

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
using CovarianceFunctions
const Kernel = CovarianceFunctions
using GaussianDistributions
using Statistics
using GaussianDistributions: optimize

# get estimate of sigma
function estimate_noise(x, y, k = Kernel.Lengthscale(Kernel.MaternP(2), .2))
    σ = 0.05 # initial guess
    σ² = optimize(σ^2, k, x, y, tol = 1e-6) # optimizing marginal likelihood
    sqrt(σ²)
end


text:
#####
##### Terms in the turbulent kinetic energy equation, all at cell centers
#####

@inline ϕ²(i, j, k, grid, ϕ) = ϕ(i, j, k, grid)^2

# Temporary way to get the vertical diffusivity for the TKE equation terms...
# Assumes that the vertical diffusivity is dominated by the CATKE contribution.
@inline shear_production(i, j, k, grid, closure, velocities, diffusivities) = zero(eltype(grid))
@inline buoyancy_flux(i, j, k, grid, closure, velocities, tracers, buoyancy, diffusivities) = zero(eltype(grid))

# Unlike the above, this fallback for dissipation is generically correct (we only want to compute dissipation once)
@inline dissipation(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs) = zero(eltype(grid))

@inline function shear_production(i, j, k, grid, closure::Union{CATKEVD, CATKEVDArray}, velocities, diffusivities)
    ∂z_u² = ℑxzᶜᵃᶜ(i, j, k, grid, ϕ², ∂zᶠᶜᶠ, velocities.u)
    ∂z_v² = ℑyzᵃᶜᶜ(i, j, k, grid, ϕ², ∂zᶜᶠᶠ, velocities.v)
    νᶻ = diffusivities.Kᵘ
    νᶻ_ijk = @inbounds νᶻ[i, j, k]
    return νᶻ_ijk * (∂z_u² + ∂z_v²)
end

@inline function buoyancy_flux(i, j, k, grid, closure::Union{CATKEVD, CATKEVDArray}, velocities, tracers, buoyancy, diffusivities)
    κᶻ = diffusivities.Kᶜ
    κᶻ_ijk = @inbounds κᶻ[i, j, k]
    N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers)
    return - κᶻ_ijk * N²
end

@inline function dissipation(i, j, k, grid, closure::Union{CATKEVD, CATKEVDArray}, velocities, tracers, buoyancy, clock, tracer_bcs)
    e = tracers.e
    FT = eltype(grid)
    @inbounds eⁱʲᵏ = e[i, j, k]
    ẽ³² = sqrt(abs(eⁱʲᵏ^3))

    ℓ = TKE_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)
    Cᴰ = closure.Cᴰ

    return Cᴰ * ẽ³² / ℓ
end

#####
##### For closure tuples...
#####

@inline shear_production(i, j, k, grid, closures::Tuple{<:Any}, velocities, diffusivities) =
    shear_production(i, j, k, grid, closures[1], velocities, diffusivities[1])

@inline shear_production(i, j, k, grid, closures::Tuple{<:Any, <:Any}, velocities, diffusivities) =
    shear_production(i, j, k, grid, closures[1], velocities, diffusivities[1]) +
    shear_production(i, j, k, grid, closures[2], velocities, diffusivities[2])

@inline shear_production(i, j, k, grid, closures::Tuple{<:Any, <:Any, <:Any}, velocities, diffusivities) =
    shear_production(i, j, k, grid, closures[1], velocities, diffusivities[1]) +
    shear_production(i, j, k, grid, closures[2], velocities, diffusivities[2]) +
    shear_production(i, j, k, grid, closures[3], velocities, diffusivities[3])

@inline buoyancy_flux(i, j, k, grid, closures::Tuple{<:Any}, velocities, tracers, buoyancy, diffusivities) =
    buoyancy_flux(i, j, k, grid, closures[1], velocities, diffusivities[1])

@inline buoyancy_flux(i, j, k, grid, closures::Tuple{<:Any, <:Any}, velocities, tracers, buoyancy, diffusivities) =
    buoyancy_flux(i, j, k, grid, closures[1], velocities, tracers, buoyancy, diffusivities[1]) +
    buoyancy_flux(i, j, k, grid, closures[2], velocities, tracers, buoyancy, diffusivities[2])

@inline buoyancy_flux(i, j, k, grid, closures::Tuple{<:Any, <:Any, <:Any}, velocities, tracers, buoyancy, diffusivities) =
    buoyancy_flux(i, j, k, grid, closures[1], velocities, tracers, buoyancy, diffusivities[1]) +
    buoyancy_flux(i, j, k, grid, closures[2], velocities, tracers, buoyancy, diffusivities[2]) +
    buoyancy_flux(i, j, k, grid, closures[3], velocities, tracers, buoyancy, diffusivities[3])

@inline dissipation(i, j, k, grid, closures::Tuple{<:Any}, velocities, tracers, buoyancy, diffusivities) =
    dissipation(i, j, k, grid, closures[1], velocities, diffusivities[1])

@inline dissipation(i, j, k, grid, closures::Tuple{<:Any, <:Any}, velocities, tracers, buoyancy, diffusivities) =
    dissipation(i, j, k, grid, closures[1], velocities, tracers, buoyancy, diffusivities[1]) +
    dissipation(i, j, k, grid, closures[2], velocities, tracers, buoyancy, diffusivities[2])

@inline dissipation(i, j, k, grid, closures::Tuple{<:Any, <:Any, <:Any}, velocities, tracers, buoyancy, diffusivities) =
    dissipation(i, j, k, grid, closures[1], velocities, tracers, buoyancy, diffusivities[1]) +
    dissipation(i, j, k, grid, closures[2], velocities, tracers, buoyancy, diffusivities[2]) +
    dissipation(i, j, k, grid, closures[3], velocities, tracers, buoyancy, diffusivities[3])


text:
const AXIOM_NODES = ["001", "000", "01"]
const DELETED = "del"

const NOISE_LEVEL_NODES = 0.6
const NOISE_LEVEL_MISSING_OUTPUT = 0.5
const EXCITATORY_MIN = 0
const EXCITATORY_MAX = 2
const INHIBITORY_THRESHOLD = 0.5

const MAX_TURN = 60

const FIGHT_THRESHOLD = 0.5
const CROSS_OVER_THRESHOLD = 0.5
const INFLUENCE_RADIUS = 2

const DISTANCE_SCALE_FACTOR = 4

const MUTATION_PROBABILITY = 0

const MAX_NEURONS = 20
const MIN_ORGANISMS = 0 # 20
const MAX_DEVELOPMENT_STEPS = 99

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# Autogenerated wrapper script for LLVM_full_jll for armv7l-linux-gnueabihf-cxx03
export clang, libclang, libclang_cpp, libllvm, liblto, llc, llvm_config, llvm_mca, mlir, opt

using Zlib_jll
JLLWrappers.@generate_wrapper_header("LLVM_full")
JLLWrappers.@declare_executable_product(clang)
JLLWrappers.@declare_library_product(libclang, "libclang.so.13")
JLLWrappers.@declare_library_product(libclang_cpp, "libclang-cpp.so.13jl")
JLLWrappers.@declare_library_product(libllvm, "libLLVM-13jl.so")
JLLWrappers.@declare_library_product(liblto, "libLTO.so.13jl")
JLLWrappers.@declare_executable_product(llc)
JLLWrappers.@declare_executable_product(llvm_config)
JLLWrappers.@declare_executable_product(llvm_mca)
JLLWrappers.@declare_library_product(mlir, "libMLIR.so.13jl")
JLLWrappers.@declare_executable_product(opt)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_executable_product(
        clang,
        "tools/clang",
    )

    JLLWrappers.@init_library_product(
        libclang,
        "lib/libclang.so",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libclang_cpp,
        "lib/libclang-cpp.so",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libllvm,
        "lib/libLLVM-13jl.so",
        nothing,
    )

    JLLWrappers.@init_library_product(
        liblto,
        "lib/libLTO.so",
        nothing,
    )

    JLLWrappers.@init_executable_product(
        llc,
        "tools/llc",
    )

    JLLWrappers.@init_executable_product(
        llvm_config,
        "tools/llvm-config",
    )

    JLLWrappers.@init_executable_product(
        llvm_mca,
        "tools/llvm-mca",
    )

    JLLWrappers.@init_library_product(
        mlir,
        "lib/libMLIR.so",
        nothing,
    )

    JLLWrappers.@init_executable_product(
        opt,
        "tools/opt",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using BRIKHEAD
using Documenter

makedocs(;
    modules=[BRIKHEAD],
    authors="Zaki A <zaki@live.ca> and contributors",
    repo="https://github.com/notZaki/BRIKHEAD.jl/blob/{commit}{path}#L{line}",
    sitename="BRIKHEAD.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://notZaki.github.io/BRIKHEAD.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/notZaki/BRIKHEAD.jl",
)


text:
#!/usr/bin/env julia

module chunk_u

#=
   Chunking Utilities - chunk_u.jl
=#

export fetch_chunk
export bound_shape, bounds, zip_bounds, intersect_bounds
export chunk_bounds, chunk_bounds2D

export vol_shape

"""

    vol_shape( bounds )

  Returns the overall shape of a subvolume which ranges
  between the two coordinates contained within the Pair bounds.
"""
function bound_shape( bounds::Pair )
  bounds.second - bounds.first + 1
end

@inline function zip_bounds( bounds::Pair )

  #need Any for possible Colon()s
  bz = Vector{Any}(length(bounds.first))

  for i in eachindex(bounds.first)
    bz[i] = bounds.first[i]:bounds.second[i]
  end

  bz
end

"""

    fetch_chunk( d, bounds::Pair, offset )

  Generalized function to fetch a 3d chunk
  from within an Nd Array, H5Dataset, etc.
"""
function fetch_chunk( d, bounds::Pair, offset=[0,0,0] )

  shifted = (bounds.first  + offset) => (bounds.second + offset);
  zipped = zip_bounds(shifted)

  while length(zipped) < length(size(d)) push!(zipped,Colon()) end

  d[zipped...]
end

"""
    chunk_bounds( vol_size, chunk_size )

  Computes the index bounds for splitting a volume
  into chunks of a maximum size. If the dimensions
  don't match to allow for perfect splitting, smaller
  are included at the end of each dimension. The offset coordinate
  is added to each resulting bound.
"""
function chunk_bounds( vol_size, chunk_size, offset=[0,0,0] )

  x_bounds = bounds1D( vol_size[1], chunk_size[1] )
  y_bounds = bounds1D( vol_size[2], chunk_size[2] )
  z_bounds = bounds1D( vol_size[3], chunk_size[3] )

  num_bounds = prod(
    (length(x_bounds), length(y_bounds), length(z_bounds))
    )
  bounds = Vector{Pair{Vector{Int},Vector{Int}}
                 }(num_bounds)
  i=0;

  for z in z_bounds
    for y in y_bounds
      for x in x_bounds
        i+=1;
        bounds[i] = [x.first,  y.first,  z.first] + offset =>
                    [x.second, y.second, z.second] + offset
      end
    end
  end

  bounds
end

"""

    chunk_bounds2D( vol_size, chunk_size, offset=[0,0] )

  Computes the index bounds for splitting a 2D image
  into chunks of a maximum size. If the dimensions
  don't match to allow for perfect splitting, smaller chunks
  are included at the end of each dimension. The offset coordinate
  is added to each resulting bound.
"""
function chunk_bounds2D( vol_size, chunk_size, offset=[0,0] )

  x_bounds = bounds1D( vol_size[1], chunk_size[1] )
  y_bounds = bounds1D( vol_size[2], chunk_size[2] )

  num_bounds = prod(
    (length(x_bounds), length(y_bounds))
  )

  bounds = Vector{Pair{Vector{Int},Vector{Int}}
                }(num_bounds)
  i=0

  for y in y_bounds
    for x in x_bounds
      i+=1;
      bounds[i] = [x.first,  y.first] + offset =>
                  [x.second, y.second] + offset
    end
  end

  bounds
end

"""

    bounds1D( full_width, step_size )

  Returns the 1D index bounds for splitting an interval
  of a given width into increments of a given step size.
  Includes smaller increments at the end if the interval
  isn't evenly divided
"""
function bounds1D( full_width, step_size )

  @assert step_size > 0
  @assert full_width > 0

  start = 1
  ending = step_size

  num_bounds = convert(Int,ceil(full_width / step_size));
  bounds = Vector{Pair{Int,Int}}( num_bounds )

  i=1;
  while ending < full_width
    bounds[i] = start => ending

    i += 1;
    start  += step_size
    ending += step_size
  end

  #last window for remainder
  bounds[end] = start => full_width

  bounds
end

"""
Tests whether a voxel coordinate is within 3d bounds
"""
function in_bounds( voxel, bounds )

  #BROADCASTING IS THE DEVIL
  # ( all(voxel .>= bounds.first) &&
  #   all(voxel .<= bounds.second)  )

  b_beg, b_end = bounds

  ( voxel[1] >= b_beg[1] &&
    voxel[2] >= b_beg[2] &&
    voxel[3] >= b_beg[3] &&

    voxel[1] <= b_end[1] &&
    voxel[2] <= b_end[2] &&
    voxel[3] <= b_end[3] )
end

function in_bounds2D( voxel, bounds )

  b_beg, b_end = bounds

  ( voxel[1] >= b_beg[1] &&
    voxel[2] >= b_beg[2] &&

    voxel[1] <= b_end[1] &&
    voxel[2] <= b_end[2] )
end

"""

    bounds( d, offset=[0,0,0] )

  Extracts the index bounds for general purpose objects
  which have size and a given offset
"""
function bounds( d, offset=[0,0,0] )
  offset + 1 => collect(size(d)[1:3]) + offset
end

"""

    intersect_bounds( bounds1, bounds2, bounds2_offset )

  Takes two index bounds defined by Pairs and returns their
  intersection. Performs no checking for validity of the result
  (i.e. the result can specify bounds with 0 or negative volume).
"""
function intersect_bounds( bounds1, bounds2, bounds2_offset )
  b1_beg, b1_end = bounds1
  b2_beg = bounds2.first  + bounds2_offset
  b2_end = bounds2.second + bounds2_offset

  max( b1_beg, b2_beg ) => min( b1_end, b2_end )
end

end#module


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

minimize a spectral function of a gram matrix of a polynomial
=#

import DynamicPolynomials

struct CentralPolyMatJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    n::Int # number of polyvars
    halfdeg::Int # half degree of random polynomials
    ext::MatSpecExt # formulation specifier
end

function build(inst::CentralPolyMatJuMP{T}) where {T <: Float64}
    (n, halfdeg) = (inst.n, inst.halfdeg)

    DynamicPolynomials.@polyvar x[1:n]
    basis = DynamicPolynomials.monomials(x, 0:halfdeg)
    L = length(basis)
    Q0 = randn(L, L)
    if is_domain_pos(inst.ext)
        # make the polynomial nonnegative
        Q0 = Q0' * Q0
    end
    Q0 .*= inv(L)
    poly = basis' * Symmetric(Q0, :U) * basis

    model = JuMP.Model()
    JuMP.@variable(model, Q_vec[1:Cones.svec_length(L)])

    # convex objective
    JuMP.@variable(model, epi)
    JuMP.@objective(model, Min, epi)
    add_homog_spectral(inst.ext, L, vcat(1.0 * epi, Q_vec), model)

    # coefficients equal
    Q = Symmetric(get_smat_U(L, 1.0 * Q_vec), :U)
    poly_eq = basis' * Q * basis - poly
    JuMP.@constraint(model, DynamicPolynomials.coefficients(poly_eq) .== 0)

    # save for use in tests
    model.ext[:Q_var] = Q

    return model
end

function test_extra(inst::CentralPolyMatJuMP{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    (stat == MOI.OPTIMAL) || return

    # check objective
    tol = eps(T)^0.2
    Q_opt = JuMP.value.(model.ext[:Q_var])
    λ = eigvals(Symmetric(Q_opt, :U))
    if is_domain_pos(inst.ext)
        @test minimum(λ) >= -tol
        λ = pos_only(λ)
    end
    obj_result = get_val(λ, inst.ext)
    @test JuMP.objective_value(model) ≈ obj_result atol=tol rtol=tol
    return
end


text:
# =========================================================================== #

# =========================================================================== #

#asked question:
#Mettre en place une heuristique de recherche locale (descente ou plus profonde descente) fondée sur deux voisinages (exemple : type “k-p exchange”).

#taking an array of lenght m
#returning another array of length m
include("zDM1_1.jl")

function calculz(x,costs,n)
    z = 0
    for i in 1:n
        z+= x[i]*costs[i]
    end

    return z
end

function select(crts,matrix,m,j)
    selec = true # juste une initialisation
    for l in 1:m
        #println(selec)
        selec = selec && (( crts[l] + matrix[l,j]) <= 1)
        #println(crts[l] + matrix[j,l] )
    end
    #println(j,selec)
    return selec
end
function kpexchange!(
    x, # Notre solution
    k::Int,  # le nombre d'objet à retirer du conteneur
    p::Int, # le nombre d'objet à rajouter dans le conteneur
    n::Int,  # taille de x
    m::Int, # nombre de contraintes
    zbest,  # solution trouvée
    z, #solution actuelle
    couts, # tableau des coûts de x
    xbest, #meilleur solution trouvée
    crts, # vecteur des contraintes (on suppose que pout tout i crts[i] =< 1)
    matrix # la matrice des contraintes
    )
    #println("===============kp exchange =====================")
    #println("x = ",x)
    sol=(xbest,zbest,x)
    if k > 0

        #println("========= k =============")
        for i in 1:n
            #println("======================")
            #println("==> élément ",i," sur ",n)
            #println("======================")
            if x[i] == 1
                #println("======================")
                #println("==> x[",i,"] = 1")
                #println("======================")
                x_prime = copy(x)
                #println("==> copie effectuée")

                x_prime[i] = 0
                #println("==> x[",i,"] passé à 0")
                for crt in 1:m
                    if matrix[crt,i] == 1
                        crts[crt] = 0
                    end
                    #println("destruction",crts)
                end
            sol = kpexchange!(x_prime,k-1,p,n,m,zbest,z-couts[i],couts, xbest,crts,matrix)
            end
        end

    elseif p > 0
        sol = kpexchange!(x,k,p-1,n,m,zbest,z,couts, xbest,crts,matrix)
        #println("========= p =============")
        for j = 1:n
            #println("======================")
            #println("==> élément ",j," sur ",n)
            #println("======================")
            if x[j] == 0
                #println("======================")
                #println("==> x[",j,"] = 0")
                #println("======================")
                # ici on regarde si on peu rajouter notre variable
                #println("crts =", crts, "" )
                x_prime = copy(sol[3])
                if select(crts,matrix,m,j)
                    #println("==> copie effectuée")
                    x_prime[j] = 1  # on rajoute la variable j à notre solution
                    #println("==> x[",j,"] passé à 1")
                    #println("======================")
                    #println("Appel récursif")
                    #z = calculz(x_prime,couts,n)
                    z = z +couts[j]
                    if z > zbest #on regarde si on a amélioré notre solution
                        xbest = x_prime
                        zbest  = z
                    else
                        z = zbest
                    end
                    sol = (zbest,xbest,x_prime)
                    #println("P CRTS :", crts)
                    for crt in 1:m # on met à jour nos contraintes
                        if matrix[crt,j] == 1
                            crts[crt] = 1
                        end
                    end
                    #println("nouvelle solution : z = ",zbest)
                end
                #println("Appel récursif p")
            end
        end

        #println(sol)
    end
    #println("zbest = ", sol[1])

    return sol
end

max_stars_repo_path: Archives/solveSPP- TOULA - ROCHET/DM1_2.jl


text:
using Test

@testset "Default Standards" begin
    @test length(getstandards(n"Rb"))>=8 # Check that the standard parse etc.
end

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
#=

  Euler #32 in Julia.

  Problem 32
  """
  We shall say that an n-digit number is pandigital if it makes use of
  all the digits 1 to n exactly once; for example, the 5-digit number,
  15234, is 1 through 5 pandigital.

  The product 7254 is unusual, as the identity, 39 × 186 = 7254,
  containing multiplicand, multiplier, and product is 1 through 9
  pandigital.

  Find the sum of all products whose multiplicand/multiplier/product
  identity can be written as a 1 through 9 pandigital.
  HINT: Some products can be obtained in more than one way so be sure
  to only include it once in your sum.
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.17556319s
function euler32a()
    prod_hash = Dict{Int,Int}()
    for a in 2:98
        as = string(a)
        if occursin("0",as) continue end
        for b in a+1:9876
            p = a*b
            if occursin("0",string(p)) continue end
            l = "$as$b$p"
            if length(l) == 9 && length(Set(l)) == 9 && !occursin("0",l)
                prod_hash[p] = 1
            end
        end
    end
    return keys(prod_hash)|>sum
end

#
function euler32b()
    prod_hash = Set()
    for a in 2:98
        # as = string(a)
        if occursin("0",string(a)) continue end
        for b in a+1:9876
            p = a*b
            if occursin("0",string(p)) continue end
            l = "$a$b$p"
            if length(l) == 9 && length(Set(l)) == 9 && !occursin("0",l)
            # if length(l) == 9 && !occursin("0",l) && length(unique(l)) == 9
                push!(prod_hash,p)
            end
        end
    end

    return prod_hash|>sum
end

function euler32c()
    prod_hash = Dict{Int,Int}()
    for a in 98:-1:2
        if a % 10 == 0 continue end
        as = string(a)
        for b in a+1:9876
            p = a*b
            if p % 10 == 0 continue end
            l = "$as$b$p"
            # l = string(a,b,p)
            if length(l) != 9 continue end
            if length(Set(l)) != 9 continue end
            if occursin("0",l) continue end
            prod_hash[p] = 1
        end
    end
    return keys(prod_hash)|>sum
end

function euler32d()
    prod_hash = Set{Int64}()
    ais = 2:98
    as = string.(ais)
    for (a,ai) in zip(as,ais)
        # if !occursin("0",a)
            for b in ai+1:9876
                p = ai*b
                if !occursin("0",string(p))
                    l = string(a,b,p)
                    if length(l) == 9 && length(Set(l)) == 9 && !occursin("0",l)
                        push!(prod_hash,p)
                    end
                end
            end
        # end
    end
    return sum(prod_hash) # keys(prod_hash)|>sum
end

function euler32e()
    prod_hash = Set{Int64}()
    for a in 2:98
        as = string(a)
        if occursin("0",as) continue end
        for b in a+1:9876
            p = a*b
            if occursin("0",string(p)) continue end
            l = string(a,b,p)
            if length(l) == 9 && length(Set(l)) == 9 && !occursin("0",l)
                push!(prod_hash,p)
            end
        end
    end
    return sum(prod_hash)
end

# 0.25s
function euler32f()
    p = [1,2,3,4,5,6,7,8,9]
    rev = reverse(p)
    ps = []
    while p != rev
        for pos in 1:4
            p1 = undigit(p[1:pos])
            p2 = undigit(p[pos+1:5])
            p3 = undigit(p[6:9])
            if p1 * p2 == p3
                push!(ps,p3)
            end
        end
        p = next_permutation(p)
    end
    return sum(unique(ps))
end

# GC.gc(true)
# GC.gc(false)
run_euler(euler32a)
# run_euler(euler32b)
# GC.gc(true)
# GC.gc(false)
# run_euler(euler32c)
# run_euler(euler32d)
# run_euler(euler32e)
# run_euler(euler32f)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# This script shows why we use the Stack datastructure to implement th e
# objectpools. Using Vector as a stack is also fast but that caused
# weird unexplainable bugs. Stack seems to work fine.
#

using StaticArrays, BenchmarkTools, DataStructures

##

module ObjectPools

using Base.Threads: nthreads, threadid

struct VectorPool{T}
    #        tid    stack    object
    arrays::Vector{Set{Vector{T}}}
    VectorPool{T}() where {T} = new( [ Set{Vector{T}}() for _=1:nthreads() ] )
end

function acquire!(pool::VectorPool{T}, len::Integer) where {T}
    tid = threadid()
    if !isempty(pool.arrays[tid]) > 0
        x = pop!(pool.arrays[tid])
        if len != length(x)
            resize!(x, len)
        end
        return x
    else
        return Vector{T}(undef, len)
    end
end

function release!(pool::VectorPool{T}, x::Vector{T}) where {T}
    tid = threadid()
    push!(pool.arrays[tid], x)
    return nothing
end

end

##

N = 1000
T = SVector{3, Float64}
vecpool = Vector{T}[]
setpool = Set{Vector{T}}()
stackpool = Stack{Vector{T}}()

for _=1:5
   push!(vecpool, rand(T, 1000))
   push!(setpool, rand(T, 1000))
   push!(stackpool, rand(T, 1000))
end

##

function runn(N, f, args...)
   for n = 1:N
      f(args...)
   end
end

function poppush(pool)
   x = pop!(pool)  # 137
   push!(pool, x)  # 199
end

runn(3, poppush, vecpool)
runn(3, poppush, setpool)
runn(3, poppush, stackpool)

##

@btime poppush($vecpool)
@btime poppush($setpool)
@btime poppush($stackpool)

##

Profile.clear()
let pool = pool
   @profile runn(100_000_000, poppush, pool)
end

Profile.print()


text:
#=
        ScalarQuadraticFunction -in- LessThan
        ScalarQuadraticFunction -in- GreaterThan
        ScalarQuadraticFunction -in- EqualTo

    TODO:
        - constraint sets
        - constraint functions
        - deleting constraints
=#
constrdict(m::LinQuadOptimizer, ::QCI{LE})  = cmap(m).q_less_than
constrdict(m::LinQuadOptimizer, ::QCI{GE})  = cmap(m).q_greater_than
constrdict(m::LinQuadOptimizer, ::QCI{EQ})  = cmap(m).q_equal_to

function MOI.addconstraint!(m::LinQuadOptimizer, func::Quad, set::S) where S <: Union{LE, GE, EQ}
    addquadraticconstraint!(m, func, set)
    m.last_constraint_reference += 1
    ref = QCI{S}(m.last_constraint_reference)
    dict = constrdict(m, ref)
    push!(m.qconstraint_primal_solution, NaN)
    push!(m.qconstraint_dual_solution, NaN)
    dict[ref] = get_number_quadratic_constraints(m)
    # dict[ref] = length(m.qconstraint_primal_solution)
    return ref
end

function addquadraticconstraint!(m::LinQuadOptimizer, func::Quad, set::S) where S<: Union{LE, GE, EQ}
    addquadraticconstraint!(m, func, backend_type(m,set), _getrhs(set))
end

function addquadraticconstraint!(m::LinQuadOptimizer, f::Quad, sense::Cchar, rhs::Float64)
    if abs(f.constant) > 0
        warn("Constant in quadratic function. Moving into set")
    end
    quadratic_columns_1 = [getcol(m, term.variable_index_1) for term in f.quadratic_terms]
    quadratic_columns_2 = [getcol(m, term.variable_index_2) for term in f.quadratic_terms]
    quadratic_coefficients = [term.coefficient for term in f.quadratic_terms]
    ri, ci, vi = reduce_duplicates!(
        quadratic_columns_1,
        quadratic_columns_2,
        quadratic_coefficients
    )
    affine_columns = [getcol(m, term.variable_index) for term in f.affine_terms]
    affine_coefficients = [term.coefficient for term in f.affine_terms]
    add_quadratic_constraint!(m,
        affine_columns,
        affine_coefficients,
        rhs - f.constant,
        sense,
        ri, ci, vi
    )
end

"""
    reduce_duplicates!(rows::Vector{T}, cols::Vector{T}, vals::Vector{S})

Given a matrix specified by row indices in `rows`, column indices in `cols` and
coefficients in `vals`, return new `rows`, `cols`, and `vals` vectors with
duplicate elements summed and any coefficients in the lower triangle moved to
the upper triangle.

This function swaps element `i` in `rows` and `cols` if `rows[i]>cols[i]`.

# Examples
```jldoctest
julia> reduce_duplicates!(
    [1,   2, 2, 2],  # rows
    [1,   1, 2, 2],  # cols
    [1, 0.5, 1, 1]   # vals
    )
([1, 1, 2], [1, 2, 2], [1.0, 0.5, 2.0])
```
"""
function reduce_duplicates!(rows::Vector{T}, cols::Vector{T}, vals::Vector{S}) where T where S
    @assert length(rows) == length(cols) == length(vals)
    for i in 1:length(rows)
        if rows[i] > cols[i]
            tmp = rows[i]
            rows[i] = cols[i]
            cols[i] = tmp
        end
    end
    findnz(sparse(rows, cols, vals))
end

MOI.candelete(m::LinQuadOptimizer, c::QCI{<: LinSets}) = MOI.isvalid(m, c)
function MOI.delete!(m::LinQuadOptimizer, c::QCI{<: LinSets})
    deleteconstraintname!(m, c)
    dict = constrdict(m, c)
    row = dict[c]
    delete_quadratic_constraints!(m, row, row)
    deleteat!(m.qconstraint_primal_solution, row)
    deleteat!(m.qconstraint_dual_solution, row)
    # shift all the other references
    shift_references_after_delete_quadratic!(m, row)
    delete!(dict, c)
end

#=
    Constraint set of Linear function
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{QCI{S}}) where S <: Union{LE, GE, EQ} = true
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintSet, c::QCI{S}) where S <: Union{LE, GE, EQ}
    rhs = get_quadratic_rhs(m, m[c])
    S(rhs)
end

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{QCI{IV}}) = false

#=
    Constraint function of Linear function
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{<:QCI{<: LinSets}}) = true
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintFunction, c::QCI{<: LinSets})
    # TODO more efficiently
    colidx, coefs, Q = get_quadratic_constraint(m, m[c])
    affine_terms = map(
        (v,c)->MOI.ScalarAffineTerm{Float64}(c, m.variable_references[v]),
        colidx,
        coefs
    )
    rows = rowvals(Q)
    vals = nonzeros(Q)
    nrows, ncols = size(Q)
    quadratic_terms = MOI.ScalarQuadraticTerm{Float64}[]
    sizehint!(quadratic_terms, length(vals))
    for i = 1:ncols
        for j in nzrange(Q, i)
            row = rows[j]
            val = vals[j]
            push!(quadratic_terms, MOI.ScalarQuadraticTerm{Float64}(2*val, m.variable_references[row], m.variable_references[i]))
        end
    end
    Quad(affine_terms, quadratic_terms, 0.0) # constant was moved into set
end


text:
#!/bin/bash
# -*- mode: julia -*-
#=
JULIA="${JULIA:-julia --color=yes --startup-file=no}"
exec ${JULIA} "$@" "${BASH_SOURCE[0]}"
=#

JL_PKG = "InteractiveCodeSearch"
VERSION >= v"0.7.0-DEV.5183" && using Pkg
Pkg.test(JL_PKG; coverage=true)


text:
#!/usr/bin/env julia

function solve(case, N, A)
    pattern = A[2] - A[1]
    longest = 2
    current = 2
    for i in 3:N
        diff = A[i] - A[i - 1]
        if diff == pattern
            current += 1
        else
            pattern = diff
            longest = max(current, longest)
            current = 2
        end
    end
    longest = max(current, longest)
    println("Case #$case: $longest")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        N = parse(Int, readline(stdin))
        A = parse.(Int, stdin |> readline |> split)
        solve(case, N, A)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
for (syevd, elty, relty) in    ((:ssyevd, :Float32, :Float32),
                            (:dsyevd, :Float64, :Float64)),
                            interface in (:Matrix, :CuMatrix)
    @eval begin
        function magma_syevd!(jobz::AbstractChar, uplo::AbstractChar, A::$interface{$elty})
            A = Matrix{$elty}(A)

            n = checksquare(A)
            lda     = max(1, stride(A, 2))
            W = similar(A, $relty, n)
            work = Vector{$elty}(undef, 1)
            lwork = -1
            iwork = Vector{Int32}(undef, 1)
            liwork = -1
            info    = Ref{Int}()

            jobz_magma      = char_to_magmaInt(jobz)
            uplo_magma     = char_to_magmaIntUplo(uplo)

            for i in 1:2
                func = eval(@magmafunc($syevd))
                func(
                    jobz_magma,
                    uplo_magma,
                    n,
                    A,
                    lda,
                    W,
                    work,
                    lwork,
                    iwork,
                    liwork,
                    info
                )

                if i == 1
                    lwork = ceil(Int, nextfloat(real(work[1])))
                    resize!(work, lwork)
                    liwork = Int(iwork[1])
                    resize!(iwork, liwork)
                end
            end

            jobz == 'V' ? (W, A) : W
        end

        function magma_syevd_m!(ngpu::Int, jobz::AbstractChar, uplo::AbstractChar, A::$interface{$elty})
            A = Matrix{$elty}(A)

            n = checksquare(A)
            lda     = max(1, stride(A, 2))
            W = similar(A, $relty, n)
            work = Vector{$elty}(undef, 1)
            lwork = -1
            iwork = Vector{Int32}(undef, 1)
            liwork = -1
            info    = Ref{Int}()

            jobz_magma      = char_to_magmaInt(jobz)
            uplo_magma     = char_to_magmaIntUplo(uplo)

            for i in 1:2
                func = eval(@magmafunc_m($syevd))
                func(
                    ngpu,
                    jobz_magma,
                    uplo_magma,
                    n,
                    A,
                    lda,
                    W,
                    work,
                    lwork,
                    iwork,
                    liwork,
                    info
                )

                if i == 1
                    lwork = ceil(Int, nextfloat(real(work[1])))
                    resize!(work, lwork)
                    liwork = Int(iwork[1])
                    resize!(iwork, liwork)
                end
            end

            jobz == 'V' ? (W, A) : W
        end
    end
end

for (syevd, elty, relty) in    ((:cheevd, :ComplexF32, :Float32),
                            (:zheevd, :ComplexF64, :Float64)),
                            interface in (:Matrix, :CuMatrix)
    @eval begin
        function magma_syevd!(jobz::AbstractChar, uplo::AbstractChar, A::$interface{$elty})
            A = Matrix{$elty}(A)

            n = checksquare(A)
            lda     = max(1, stride(A, 2))
            W = similar(A, $relty, n)
            work = Vector{$elty}(undef, 1)
            lwork = -1
            rwork = Vector{$relty}(undef, 1)
            lrwork = -1
            iwork = Vector{Int32}(undef, 1)
            liwork = -1
            info    = Ref{Int}()

            jobz_magma      = char_to_magmaInt(jobz)
            uplo_magma     = char_to_magmaIntUplo(uplo)

            for i in 1:2
                func = eval(@magmafunc($syevd))
                func(
                    jobz_magma,
                    uplo_magma,
                    n,
                    A,
                    lda,
                    W,
                    work,
                    lwork,
                    rwork,
                    lrwork,
                    iwork,
                    liwork,
                    info
                )

                if i == 1
                    lwork = ceil(Int, nextfloat(real(work[1])))
                    resize!(work, lwork)
                    lrwork = ceil(Int, nextfloat(real(rwork[1])))
                    resize!(rwork, lrwork)
                    liwork = Int(iwork[1])
                    resize!(iwork, liwork)
                end
            end

            jobz == 'V' ? (W, A) : W
        end

        function magma_syevd_m!(ngpu::Int, jobz::AbstractChar, uplo::AbstractChar, A::$interface{$elty})
            A = Matrix{$elty}(A)

            n = checksquare(A)
            lda     = max(1, stride(A, 2))
            W = similar(A, $relty, n)
            work = Vector{$elty}(undef, 1)
            lwork = -1
            rwork = Vector{$relty}(undef, 1)
            lrwork = -1
            iwork = Vector{Int32}(undef, 1)
            liwork = -1
            info    = Ref{Int}()

            jobz_magma      = char_to_magmaInt(jobz)
            uplo_magma     = char_to_magmaIntUplo(uplo)

            for i in 1:2
                func = eval(@magmafunc_m($syevd))
                func(
                    ngpu,
                    jobz_magma,
                    uplo_magma,
                    n,
                    A,
                    lda,
                    W,
                    work,
                    lwork,
                    rwork,
                    lrwork,
                    iwork,
                    liwork,
                    info
                )

                if i == 1
                    lwork = ceil(Int, nextfloat(real(work[1])))
                    resize!(work, lwork)
                    lrwork = ceil(Int, nextfloat(real(rwork[1])))
                    resize!(rwork, lrwork)
                    liwork = Int(iwork[1])
                    resize!(iwork, liwork)
                end
            end

            jobz == 'V' ? (W, A) : W
        end
    end
end


text:
@testset "shortest path" begin
    using STMO
    using STMO.ShortestPath

    graph = Dict('A' => [(0.5, 'B'), (1.5, 'H')],
        'B' => [(1.0, 'C')],
        'C' => [(1.0, 'D'), (1.0, 'E')],
        'D' => [(2.0, 'F')],
        'E' => [(1.0, 'G')],
        'F' => [(0.1, 'E')],
        'G' => [(1.0, 'F')],
        'H' => [(0.9, 'E')])

    distances, previous = dijkstra(graph, 'A')

    @test distances['D'] ≈ 2.5
    @test previous['E'] == 'H'

    path = reconstruct_path(previous, 'A', 'H')
    @test first(path) == 'A'
    @test last(path) == 'H'

    path_d, cost_AH = dijkstra(graph, 'A', 'H')
    @test path_d == path
    @test cost_AH ≈ distances['H']

    path_as, cost_as = a_star(graph, 'A', 'H', x -> 0.0)
    @test path_as == path
    @test cost_as ≈ distances['H']

end