text:
export TDLearner, DifferentialTDLearner, TDλReturnLearner

using LinearAlgebra:dot
using Distributions:pdf

"""
    TDLearner(approximator::Tapp, γ::Float64, optimizer::Float64; n::Int=0) where {Tapp<:AbstractVApproximator}
    TDLearner(approximator::Tapp, γ::Float64, optimizer::Float64; n::Int=0, method::Symbol=:SARSA) where {Tapp<:AbstractApproximator}

The `TDLearner`(Temporal Difference Learner) use the latest `n` step experiences to update the `approximator`.

# Args

- `approximator`::[`AbstractApproximator`](@ref)
- Note that `n` starts with `0`, which means looking forward for the next `n` steps.
- `γ` is the discount rate of experience.
- `optimizer` is the learning rate.
- `method`. For [`VApproximator`](@ref), the only supported update method is `:SRS`, which means only **S**tates, **R**ewards and next_**S**ates are used to update the `approximator`. For [`QApproximator`](@ref), the following methods are supported:
    - `:SARS` (aka Q-Learning)
    - `:SARSA`
    - `:ExpectedSARSA`
"""
mutable struct TDLearner{Tapp<:AbstractApproximator,method,O} <: AbstractLearner
    approximator::Tapp
    γ::Float64
    optimizer::O
    n::Int

    function TDLearner(
        ;
        approximator::Tapp,
        γ = 1.0,
        optimizer = Descent(1.0),
        n = 0,
        method::Symbol = :SARSA,
    ) where {Tapp<:AbstractApproximator}
        if ApproximatorStyle(approximator) === VApproximator()
            if method != :SRS
                throw(ArgumentError("method [$method] is unsupported for a value approximator"))
            end
        elseif ApproximatorStyle(approximator) === QApproximator()
            if !in(method, [:SARSA, :SARS, :ExpectedSARSA])
                throw(ArgumentError("Supported methods are $supported_methods , your input is $method"))
            end
        else
            throw(ArgumentError("unknown approximator style"))
        end
        new{Tapp,method,typeof(optimizer)}(approximator, γ, optimizer, n)
    end
end

(learner::TDLearner)(obs) = learner.approximator(get_state(obs))
(learner::TDLearner)(obs, a) = learner.approximator(get_state(obs))[a]

RLBase.update!(learner::TDLearner{T, M}, experience::NamedTuple) where {T, M} = update!(learner, ApproximatorStyle(learner.approximator), Val(M), experience)

extract_experience(t::AbstractTrajectory, learner::TDLearner{T, M}) where {T, M} = extract_experience(t, learner, ApproximatorStyle(learner.approximator), Val(M))

#####
# SARSA
#####

function RLBase.update!(
    learner::TDLearner,
    ::QApproximator,
    ::Val{:SARSA},
    transitions::NamedTuple{(:states, :actions, :rewards, :terminals, :next_states, :next_actions)}
)
    states, actions, rewards, terminals, next_states, next_actions = transitions
    n, γ, Q, optimizer = learner.n, learner.γ, learner.approximator, learner.optimizer

    if length(terminals) > 0 && terminals[end]
        @views gains = discount_rewards(rewards[max(end - n, 1):end], γ)  # n starts with 0
        for (i, G) in enumerate(gains)
            @views s, a = states[end-length(gains)+i], actions[end-length(gains)+i]
            update!(Q, (s, a) => apply!(optimizer, (s, a), G - Q(s, a)))
        end
    else
        if length(states) ≥ (n + 1)  # n starts with 0
            @views s, a, s′, a′ = states[end-n],
                actions[end-n],
                next_states[end],
                next_actions[end]
            @views G = discount_rewards_reduced(rewards[end-n:end], γ) +
                       γ^(n + 1) * Q(s′, a′)
            update!(Q, (s, a) => apply!(optimizer, (s, a), G - Q(s, a)))
        end
    end
end

function extract_experience(
    t::AbstractTrajectory,
    learner::TDLearner,
    ::QApproximator,
    ::Val{:SARSA},
)
    n, N = learner.n, length(t)
    # !!! n starts with 0
    if N > 0
        (
            states = select_last_dim(get_trace(t, :state), max(1, N-n):N),
            actions = select_last_dim(get_trace(t, :action), max(1, N-n):N),
            rewards = select_last_dim(get_trace(t, :reward), max(1, N-n):N),
            terminals = select_last_dim(get_trace(t, :terminal), max(1, N-n):N),
            next_states = select_last_dim(get_trace(t, :next_state), max(1, N-n):N),
            next_actions = select_last_dim(get_trace(t, :next_action), max(1, N-n):N),
        )
    else
        nothing
    end
end

#####
# ExpectedSARSA
#####

function RLBase.update!(
    learner::TDLearner,
    ::QApproximator,
    ::Val{:ExpectedSARSA},
    transitions::NamedTuple{(:states, :actions, :rewards, :terminals, :next_states, :prob_of_next_actions)}
)
    states, actions, rewards, terminals, next_states, probs_of_a′ = transitions
    n, γ, Q, optimizer = learner.n, learner.γ, learner.approximator, learner.optimizer

    if length(terminals) > 0 && terminals[end]
        @views gains = discount_rewards(rewards[max(end - n, 1):end], γ)  # n starts with 0
        for (i, G) in enumerate(gains)
            @views s, a = states[end-length(gains)+i], actions[end-length(gains)+i]
            update!(Q, (s, a) => apply!(optimizer, (s, a), G - Q(s, a)))
        end
    else
        if length(states) ≥ (n + 1)  # n starts with 0
            @views s, a, s′ = states[end-n], actions[end-n], next_states[end]
            @views G = discount_rewards_reduced(rewards[end-n:end], γ) +
                       γ^(n + 1) * dot(Q(s′), probs_of_a′)
            update!(Q, (s, a) => apply!(optimizer, (s, a), G - Q(s, a)))
        end
    end
end

function extract_experience(
    t::AbstractTrajectory,
    policy::QBasedPolicy{<:TDLearner{<:AbstractApproximator,:ExpectedSARSA}}
)
    if length(t) > 0
        # !!! n starts with 0
        n, N = policy.learner.n, length(t)
        (
            states = select_last_dim(get_trace(t, :state), max(1, N-n):N),
            actions = select_last_dim(get_trace(t, :action), max(1, N-n):N),
            rewards = select_last_dim(get_trace(t, :reward), max(1, N-n):N),
            terminals = select_last_dim(get_trace(t, :terminal), max(1, N-n):N),
            next_states = select_last_dim(get_trace(t, :next_state), max(1, N-n):N),
            prob_of_next_actions = pdf(get_prob(policy, (state=select_last_frame(get_trace(t, :next_state)),)))
        )
    else
        nothing
    end
end

#####
# SARS
#####

function RLBase.update!(
    learner::TDLearner,
    ::QApproximator,
    ::Val{:SARS},
    transitions::NamedTuple{(:states, :actions, :rewards, :terminals, :next_states)}
)
    states, actions, rewards, terminals, next_states = transitions
    n, γ, Q, optimizer = learner.n, learner.γ, learner.approximator, learner.optimizer

    if length(terminals) > 0 && terminals[end]
        @views gains = discount_rewards(rewards[max(end - n, 1):end], γ)  # n starts with 0
        for (i, G) in enumerate(gains)
            @views s, a = states[end-length(gains)+i], actions[end-length(gains)+i]
            update!(Q, (s, a) => apply!(optimizer, (s, a), G - Q(s, a)))
        end
    else
        if length(states) ≥ (n + 1)  # n starts with 0
            @views s, a, s′ = states[end-n], actions[end-n], next_states[end]
            @views G = discount_rewards_reduced(rewards[end-n:end], γ) +
                       γ^(n + 1) * maximum(Q(s′))  # n starts with 0
            update!(Q, (s, a) => apply!(optimizer, (s, a), G - Q(s, a)))
        end
    end
end

function extract_experience(
    t::AbstractTrajectory,
    learner::TDLearner,
    ::QApproximator,
    ::Val{:SARS},
)
    n, N = learner.n, length(t)
    if length(t) > 0
        (
            states = select_last_dim(get_trace(t, :state), max(1, N-n):N),
            actions = select_last_dim(get_trace(t, :action), max(1, N-n):N),
            rewards = select_last_dim(get_trace(t, :reward), max(1, N-n):N),
            terminals = select_last_dim(get_trace(t, :terminal), max(1, N-n):N),
            next_states = select_last_dim(get_trace(t, :next_state), max(1, N-n):N),
        )
    else
        nothing
    end
end

function RLBase.update!(
    learner::TDLearner{<:AbstractApproximator,:SARS},
    model::Union{TimeBasedSampleModel,ExperienceBasedSampleModel},
    t::AbstractTrajectory,
    plan_step::Int,
)
    @assert learner.n == 0 "n must be 0 here"
    for _ = 1:plan_step
        transitions = extract_experience(model, learner)
        if !isnothing(transitions)
            update!(learner, transitions)
        end
    end
end

function extract_experience(
    model::Union{ExperienceBasedSampleModel,TimeBasedSampleModel},
    learner::TDLearner{<:AbstractApproximator,:SARS},
)
    if length(model.experiences) > 0
        s = sample(model)
        (
            states = [s[1]],
            actions = [s[2]],
            rewards = [s[3]],
            terminals = [s[4]],
            next_states = [s[5]],
        )
    else
        nothing
    end
end

function RLBase.get_priority(learner::TDLearner{<:AbstractApproximator,:SARS}, transition::Tuple)
    s, a, r, d, s′ = transition
    γ, Q, opt = learner.γ, learner.approximator, learner.optimizer
    error = d ? apply!(opt, (s, a), r - Q(s, a)) :
            apply!(opt, (s, a), r + γ^(learner.n + 1) * maximum(Q(s′)) - Q(s, a))
    abs(error)
end

function RLBase.update!(
    learner::TDLearner{<:AbstractApproximator,:SARS},
    model::PrioritizedSweepingSampleModel,
    t::AbstractTrajectory,
    plan_step::Int,
)
    for _ = 1:plan_step
        # @assert learner.n == 0 "n must be 0 here"
        transitions = extract_experience(model, learner)
        if !isnothing(transitions)
            update!(learner, transitions)
            s, _, _, _, _ = transitions
            s = s[]  # length(s) is assumed to be 1
            for (s̄, ā, r̄, d̄) in model.predecessors[s]
                P = get_priority(learner, (s̄, ā, r̄, d̄, s))
                if P ≥ model.θ
                    model.PQueue[(s̄, ā)] = P
                end
            end
        end
    end
end

function extract_experience(
    model::PrioritizedSweepingSampleModel,
    learner::TDLearner{<:AbstractApproximator,:SARS},
)
    if length(model.PQueue) > 0
        s = sample(model)
        (
            states = [s[1]],
            actions = [s[2]],
            rewards = [s[3]],
            terminals = [s[4]],
            next_states = [s[5]],
        )
    else
        nothing
    end
end

#####
# SARS DoubleLearner
#####

function RLBase.update!(
    learner::DoubleLearner{T, T},
    transitions::NamedTuple{(:states, :actions, :rewards, :terminals, :next_states)},
) where {T<:TDLearner{<:AbstractApproximator,:SARS}}

    if rand(learner.rng, Bool)
        learner, target_learner = learner.L1, learner.L2
    else
        learner, target_learner = learner.L2, learner.L1
    end

    states, actions, rewards, terminals, next_states = transitions
    n, γ, Q, Qₜ, optimizer = learner.n,
        learner.γ,
        learner.approximator,
        target_learner.approximator,
        learner.optimizer

    if length(terminals) > 0 && terminals[end]
        @views gains = discount_rewards(rewards[max(end - n, 1):end], γ)  # n starts with 0
        for (i, G) in enumerate(gains)
            @views s, a = states[end-length(gains)+i], actions[end-length(gains)+i]
            update!(Q, (s, a) => apply!(optimizer, (s, a), G - Q(s, a)))
        end
    else
        if length(states) ≥ (n + 1)  # n starts with 0
            @views s, a, s′ = states[end-n], actions[end-n], next_states[end]
            @views G = discount_rewards_reduced(rewards[end-n:end], γ) +
                       γ^(n + 1) * Qₜ(s′, argmax(Q(s′)))
            update!(Q, (s, a) => apply!(optimizer, (s, a), G - Q(s, a)))
        end
    end
end

#####
# SRS
#####

RLBase.update!(
    learner::TDLearner,
    ::VApproximator,
    ::Val{:SRS},
    transitions::NamedTuple{(:states, :rewards, :terminals, :next_states)}
) = update!(
    learner,
    (
        states = transitions.states,
        actions = nothing,
        rewards = transitions.rewards,
        terminals = transitions.terminals,
        next_states = transitions.next_states,
        weights = nothing
    )
)

function RLBase.update!(
    learner::TDLearner,
    ::VApproximator,
    ::Val{:SRS},
    transitions::NamedTuple{(:states, :actions, :rewards, :terminals, :next_states, :weights)}
)
    states, actions, rewards, terminals, next_states, weights = transitions
    n, γ, V, optimizer = learner.n, learner.γ, learner.approximator, learner.optimizer

    if length(terminals) > 0 && terminals[end]
        @views gains = discount_rewards(rewards[max(end - n, 1):end], γ)  # n starts with 0
        cum_weights = isnothing(weights) ? nothing : cumprod!(reverse(weights))
        for (i, G) in enumerate(gains)
            @views s = states[end-length(gains)+i]
            if isnothing(cum_weights)
                update!(V, s => apply!(optimizer, s, G - V(s)))
            else
                update!(V, s => apply!(optimizer, s, cum_weights[i] * (G - V(s))))
            end
        end
    else
        if length(states) ≥ (n + 1)  # n starts with 0
            @views G = discount_rewards_reduced(rewards[end-n:end], γ) +
                       γ^(n + 1) * V(next_states[end])
            @views s = states[end-n]
            w = isnothing(weights) ? 1.0 : reduce(*, weights)
            update!(V, s => apply!(optimizer, s, w * (G - V(s))))
        end
    end
end

function extract_experience(
    t::AbstractTrajectory,
    learner::TDLearner,
    ::VApproximator,
    ::Val{:SRS},
)
    n, N = learner.n, length(t)
    if length(t) > 0
        (
            states = select_last_dim(get_trace(t, :state), max(1, N-n):N),
            rewards = select_last_dim(get_trace(t, :reward), max(1, N-n):N),
            terminals = select_last_dim(get_trace(t, :terminal), max(1, N-n):N),
            next_states = select_last_dim(get_trace(t, :next_state), max(1, N-n):N),
        )
    else
        nothing
    end
end

function extract_experience(
    t::AbstractTrajectory,
    π::OffPolicy{<:VBasedPolicy{<:TDLearner{<:AbstractApproximator,:SRS}}},
)
    transitions = extract_experience(t, π.π_target.learner)
    if isnothing(transitions)
        nothing
    else
        n, N = π.π_target.learner.n, length(t)
        (
            states = transitions.states,
            actions = select_last_dim(get_trace(t, :action), max(1, N-n):N),
            rewards = transitions.rewards,
            terminals = transitions.terminals,
            next_states = transitions.next_states,
        )
    end
end

#####
# DifferentialTDLearner
#####

"""
    DifferentialTDLearner(;approximator::A, α::Float64, β::Float64, R̄::Float64 = 0.0, n::Int = 0)
"""
Base.@kwdef mutable struct DifferentialTDLearner{A<:AbstractApproximator} <: AbstractLearner
    approximator::A
    α::Float64
    β::Float64
    R̄::Float64 = 0.0
    n::Int = 0
end

(learner::DifferentialTDLearner)(obs) = learner.approximator(get_state(obs))
(learner::DifferentialTDLearner)(obs, a) = learner.approximator(get_state(obs), a)

function RLBase.update!(learner::DifferentialTDLearner, transition::NamedTuple)
    states, actions, rewards, terminals, next_states, next_actions = transition
    n, α, β, Q = learner.n, learner.α, learner.β, learner.approximator
    if length(states) ≥ n + 1
        s, a = states[1], actions[1]
        s′, a′ = next_states[end], next_actions[end]
        δ = sum(r -> r - learner.R̄, rewards) + Q(s′, a′) - Q(s, a)
        learner.R̄ += β * δ
        update!(Q, (s, a) => α * δ)
    end
end

function extract_experience(
    t::AbstractTrajectory,
    learner::DifferentialTDLearner,
)
    n, N = learner.n, length(t)
    # !!! n starts with 0
    if N > 0
        (
            states = select_last_dim(get_trace(t, :state), max(1, N-n):N),
            actions = select_last_dim(get_trace(t, :action), max(1, N-n):N),
            rewards = select_last_dim(get_trace(t, :reward), max(1, N-n):N),
            terminals = select_last_dim(get_trace(t, :terminal), max(1, N-n):N),
            next_states = select_last_dim(get_trace(t, :next_state), max(1, N-n):N),
            next_actions = select_last_dim(get_trace(t, :next_action), max(1, N-n):N),
        )
    else
        nothing
    end
end

#####
# TDλReturnLearner
#####
"""
    TDλReturnLearner(;approximator::Tapp, γ::Float64 = 1.0, α::Float64, λ::Float64)
"""
Base.@kwdef struct TDλReturnLearner{Tapp<:AbstractApproximator} <: AbstractLearner
    approximator::Tapp
    γ::Float64 = 1.0
    α::Float64
    λ::Float64
end

(learner::TDλReturnLearner)(obs) = learner.approximator(obs)
(learner::TDλReturnLearner)(obs, a) = learner.approximator(obs, a)

function extract_experience(
    t::AbstractTrajectory,
    learner::TDλReturnLearner{<:AbstractApproximator},
)
    if length(t) > 0 && get_trace(t, :terminal)[end]
        (
            states = get_trace(t, :state),
            rewards = get_trace(t, :reward),
            terminals = get_trace(t, :terminal),
            next_states = get_trace(t, :next_state),
        )
    else
        nothing
    end
end

function RLBase.update!(learner::TDλReturnLearner, transition::NamedTuple)
    λ, γ, V, α = learner.λ, learner.γ, learner.approximator, learner.α
    states, rewards, terminals, next_states = transition
    T = length(states)
    for t = 1:T
        G = 0.0
        for n = 1:(T-t)
            G += λ^(n - 1) *
                 (discount_rewards_reduced(@view(rewards[t:t+n-1]), γ) +
                  γ^n * V(next_states[t+n-1]))
        end
        G *= 1 - λ
        G += λ^(T - t) *
             (discount_rewards_reduced(@view(rewards[t:T]), γ) +
              γ^(T - t + 1) * V(next_states[T]))
        sₜ = states[t]
        update!(V, sₜ => α * (G - V(sₜ)))
    end
end

ext: jl


text:
global const PetscErrorType = Int32
global const PETSC_ERROR_INITIAL = Int32( 0 )
global const PETSC_ERROR_REPEAT = Int32( 1 )
global const PETSC_ERROR_IN_CXX = Int32( 2 )

global const KSPConvergedReason = Cint
global const KSP_CONVERGED_RTOL_NORMAL = Cint( 1 )
global const KSP_CONVERGED_ATOL_NORMAL = Cint( 9 )
global const KSP_CONVERGED_RTOL = Cint( 2 )
global const KSP_CONVERGED_ATOL = Cint( 3 )
global const KSP_CONVERGED_ITS = Cint( 4 )
global const KSP_CONVERGED_CG_NEG_CURVE = Cint( 5 )
global const KSP_CONVERGED_CG_CONSTRAINED = Cint( 6 )
global const KSP_CONVERGED_STEP_LENGTH = Cint( 7 )
global const KSP_CONVERGED_HAPPY_BREAKDOWN = Cint( 8 )
global const KSP_DIVERGED_NULL = Cint( -2 )
global const KSP_DIVERGED_ITS = Cint( -3 )
global const KSP_DIVERGED_DTOL = Cint( -4 )
global const KSP_DIVERGED_BREAKDOWN = Cint( -5 )
global const KSP_DIVERGED_BREAKDOWN_BICG = Cint( -6 )
global const KSP_DIVERGED_NONSYMMETRIC = Cint( -7 )
global const KSP_DIVERGED_INDEFINITE_PC = Cint( -8 )
global const KSP_DIVERGED_NANORINF = Cint( -9 )
global const KSP_DIVERGED_INDEFINITE_MAT = Cint( -10 )
global const KSP_DIVERGED_PCSETUP_FAILED = Cint( -11 )
global const KSP_CONVERGED_ITERATING = Cint( 0 )

global const MatOption = Cint
global const MAT_OPTION_MIN = Cint( -3 )
global const MAT_NEW_NONZERO_LOCATION_ERR = Cint( 11 )
global const MAT_UNUSED_NONZERO_LOCATION_ERR = Cint( -2 )
global const MAT_NEW_NONZERO_ALLOCATION_ERR = Cint( 19 )
global const MAT_ROW_ORIENTED = Cint( -1 )
global const MAT_SYMMETRIC = Cint( 1 )
global const MAT_STRUCTURALLY_SYMMETRIC = Cint( 2 )
global const MAT_NEW_DIAGONALS = Cint( 3 )
global const MAT_IGNORE_OFF_PROC_ENTRIES = Cint( 4 )
global const MAT_USE_HASH_TABLE = Cint( 5 )
global const MAT_KEEP_NONZERO_PATTERN = Cint( 6 )
global const MAT_IGNORE_ZERO_ENTRIES = Cint( 7 )
global const MAT_USE_INODES = Cint( 8 )
global const MAT_HERMITIAN = Cint( 9 )
global const MAT_SYMMETRY_ETERNAL = Cint( 10 )
global const MAT_IGNORE_LOWER_TRIANGULAR = Cint( 12 )
global const MAT_ERROR_LOWER_TRIANGULAR = Cint( 13 )
global const MAT_GETROW_UPPERTRIANGULAR = Cint( 14 )
global const MAT_SPD = Cint( 15 )
global const MAT_NO_OFF_PROC_ZERO_ROWS = Cint( 16 )
global const MAT_NO_OFF_PROC_ENTRIES = Cint( 17 )
global const MAT_NEW_NONZERO_LOCATIONS = Cint( 18 )
global const MAT_OPTION_MAX = Cint( 21 )

global const MatInfoType = Int32
global const MAT_LOCAL = Int32( 1 )
global const MAT_GLOBAL_MAX = Int32( 2 )
global const MAT_GLOBAL_SUM = Int32( 3 )

global const KSPNormType = Cint
global const KSP_NORM_DEFAULT = Cint( -1 )
global const KSP_NORM_NONE = Cint( 0 )
global const KSP_NORM_PRECONDITIONED = Cint( 1 )
global const KSP_NORM_UNPRECONDITIONED = Cint( 2 )
global const KSP_NORM_NATURAL = Cint( 3 )

global const MatOperation = UInt32
global const MATOP_SET_VALUES = UInt32( 0 )
global const MATOP_GET_ROW = UInt32( 1 )
global const MATOP_RESTORE_ROW = UInt32( 2 )
global const MATOP_MULT = UInt32( 3 )
global const MATOP_MULT_ADD = UInt32( 4 )
global const MATOP_MULT_TRANSPOSE = UInt32( 5 )
global const MATOP_MULT_TRANSPOSE_ADD = UInt32( 6 )
global const MATOP_SOLVE = UInt32( 7 )
global const MATOP_SOLVE_ADD = UInt32( 8 )
global const MATOP_SOLVE_TRANSPOSE = UInt32( 9 )
global const MATOP_SOLVE_TRANSPOSE_ADD = UInt32( 10 )
global const MATOP_LUFACTOR = UInt32( 11 )
global const MATOP_CHOLESKYFACTOR = UInt32( 12 )
global const MATOP_SOR = UInt32( 13 )
global const MATOP_TRANSPOSE = UInt32( 14 )
global const MATOP_GETINFO = UInt32( 15 )
global const MATOP_EQUAL = UInt32( 16 )
global const MATOP_GET_DIAGONAL = UInt32( 17 )
global const MATOP_DIAGONAL_SCALE = UInt32( 18 )
global const MATOP_NORM = UInt32( 19 )
global const MATOP_ASSEMBLY_BEGIN = UInt32( 20 )
global const MATOP_ASSEMBLY_END = UInt32( 21 )
global const MATOP_SET_OPTION = UInt32( 22 )
global const MATOP_ZERO_ENTRIES = UInt32( 23 )
global const MATOP_ZERO_ROWS = UInt32( 24 )
global const MATOP_LUFACTOR_SYMBOLIC = UInt32( 25 )
global const MATOP_LUFACTOR_NUMERIC = UInt32( 26 )
global const MATOP_CHOLESKY_FACTOR_SYMBOLIC = UInt32( 27 )
global const MATOP_CHOLESKY_FACTOR_NUMERIC = UInt32( 28 )
global const MATOP_SETUP_PREALLOCATION = UInt32( 29 )
global const MATOP_ILUFACTOR_SYMBOLIC = UInt32( 30 )
global const MATOP_ICCFACTOR_SYMBOLIC = UInt32( 31 )
global const MATOP_DUPLICATE = UInt32( 34 )
global const MATOP_FORWARD_SOLVE = UInt32( 35 )
global const MATOP_BACKWARD_SOLVE = UInt32( 36 )
global const MATOP_ILUFACTOR = UInt32( 37 )
global const MATOP_ICCFACTOR = UInt32( 38 )
global const MATOP_AXPY = UInt32( 39 )
global const MATOP_GET_SUBMATRICES = UInt32( 40 )
global const MATOP_INCREASE_OVERLAP = UInt32( 41 )
global const MATOP_GET_VALUES = UInt32( 42 )
global const MATOP_COPY = UInt32( 43 )
global const MATOP_GET_ROW_MAX = UInt32( 44 )
global const MATOP_SCALE = UInt32( 45 )
global const MATOP_SHIFT = UInt32( 46 )
global const MATOP_DIAGONAL_SET = UInt32( 47 )
global const MATOP_ZERO_ROWS_COLUMNS = UInt32( 48 )
global const MATOP_SET_RANDOM = UInt32( 49 )
global const MATOP_GET_ROW_IJ = UInt32( 50 )
global const MATOP_RESTORE_ROW_IJ = UInt32( 51 )
global const MATOP_GET_COLUMN_IJ = UInt32( 52 )
global const MATOP_RESTORE_COLUMN_IJ = UInt32( 53 )
global const MATOP_FDCOLORING_CREATE = UInt32( 54 )
global const MATOP_COLORING_PATCH = UInt32( 55 )
global const MATOP_SET_UNFACTORED = UInt32( 56 )
global const MATOP_PERMUTE = UInt32( 57 )
global const MATOP_SET_VALUES_BLOCKED = UInt32( 58 )
global const MATOP_GET_SUBMATRIX = UInt32( 59 )
global const MATOP_DESTROY = UInt32( 60 )
global const MATOP_VIEW = UInt32( 61 )
global const MATOP_CONVERT_FROM = UInt32( 62 )
global const MATOP_MATMAT_MULT = UInt32( 63 )
global const MATOP_MATMAT_MULT_SYMBOLIC = UInt32( 64 )
global const MATOP_MATMAT_MULT_NUMERIC = UInt32( 65 )
global const MATOP_SET_LOCAL_TO_GLOBAL_MAP = UInt32( 66 )
global const MATOP_SET_VALUES_LOCAL = UInt32( 67 )
global const MATOP_ZERO_ROWS_LOCAL = UInt32( 68 )
global const MATOP_GET_ROW_MAX_ABS = UInt32( 69 )
global const MATOP_GET_ROW_MIN_ABS = UInt32( 70 )
global const MATOP_CONVERT = UInt32( 71 )
global const MATOP_SET_COLORING = UInt32( 72 )
global const MATOP_SET_VALUES_ADIFOR = UInt32( 74 )
global const MATOP_FD_COLORING_APPLY = UInt32( 75 )
global const MATOP_SET_FROM_OPTIONS = UInt32( 76 )
global const MATOP_MULT_CONSTRAINED = UInt32( 77 )
global const MATOP_MULT_TRANSPOSE_CONSTRAIN = UInt32( 78 )
global const MATOP_FIND_ZERO_DIAGONALS = UInt32( 79 )
global const MATOP_MULT_MULTIPLE = UInt32( 80 )
global const MATOP_SOLVE_MULTIPLE = UInt32( 81 )
global const MATOP_GET_INERTIA = UInt32( 82 )
global const MATOP_LOAD = UInt32( 83 )
global const MATOP_IS_SYMMETRIC = UInt32( 84 )
global const MATOP_IS_HERMITIAN = UInt32( 85 )
global const MATOP_IS_STRUCTURALLY_SYMMETRIC = UInt32( 86 )
global const MATOP_SET_VALUES_BLOCKEDLOCAL = UInt32( 87 )
global const MATOP_GET_VECS = UInt32( 88 )
global const MATOP_MAT_MULT = UInt32( 89 )
global const MATOP_MAT_MULT_SYMBOLIC = UInt32( 90 )
global const MATOP_MAT_MULT_NUMERIC = UInt32( 91 )
global const MATOP_PTAP = UInt32( 92 )
global const MATOP_PTAP_SYMBOLIC = UInt32( 93 )
global const MATOP_PTAP_NUMERIC = UInt32( 94 )
global const MATOP_MAT_TRANSPOSE_MULT = UInt32( 95 )
global const MATOP_MAT_TRANSPOSE_MULT_SYMBO = UInt32( 96 )
global const MATOP_MAT_TRANSPOSE_MULT_NUMER = UInt32( 97 )
global const MATOP_CONJUGATE = UInt32( 102 )
global const MATOP_SET_VALUES_ROW = UInt32( 104 )
global const MATOP_REAL_PART = UInt32( 105 )
global const MATOP_IMAGINARY_PART = UInt32( 106 )
global const MATOP_GET_ROW_UPPER_TRIANGULAR = UInt32( 107 )
global const MATOP_RESTORE_ROW_UPPER_TRIANG = UInt32( 108 )
global const MATOP_MAT_SOLVE = UInt32( 109 )
global const MATOP_GET_REDUNDANT_MATRIX = UInt32( 110 )
global const MATOP_GET_ROW_MIN = UInt32( 111 )
global const MATOP_GET_COLUMN_VECTOR = UInt32( 112 )
global const MATOP_MISSING_DIAGONAL = UInt32( 113 )
global const MATOP_GET_SEQ_NONZERO_STRUCTUR = UInt32( 114 )
global const MATOP_CREATE = UInt32( 115 )
global const MATOP_GET_GHOSTS = UInt32( 116 )
global const MATOP_GET_LOCAL_SUB_MATRIX = UInt32( 117 )
global const MATOP_RESTORE_LOCALSUB_MATRIX = UInt32( 118 )
global const MATOP_MULT_DIAGONAL_BLOCK = UInt32( 119 )
global const MATOP_HERMITIAN_TRANSPOSE = UInt32( 120 )
global const MATOP_MULT_HERMITIAN_TRANSPOSE = UInt32( 121 )
global const MATOP_MULT_HERMITIAN_TRANS_ADD = UInt32( 122 )
global const MATOP_GET_MULTI_PROC_BLOCK = UInt32( 123 )
global const MATOP_FIND_NONZERO_ROWS = UInt32( 124 )
global const MATOP_GET_COLUMN_NORMS = UInt32( 125 )
global const MATOP_INVERT_BLOCK_DIAGONAL = UInt32( 126 )
global const MATOP_GET_SUB_MATRICES_PARALLE = UInt32( 128 )
global const MATOP_SET_VALUES_BATCH = UInt32( 129 )
global const MATOP_TRANSPOSE_MAT_MULT = UInt32( 130 )
global const MATOP_TRANSPOSE_MAT_MULT_SYMBO = UInt32( 131 )
global const MATOP_TRANSPOSE_MAT_MULT_NUMER = UInt32( 132 )
global const MATOP_TRANSPOSE_COLORING_CREAT = UInt32( 133 )
global const MATOP_TRANS_COLORING_APPLY_SPT = UInt32( 134 )
global const MATOP_TRANS_COLORING_APPLY_DEN = UInt32( 135 )
global const MATOP_RART = UInt32( 136 )
global const MATOP_RART_SYMBOLIC = UInt32( 137 )
global const MATOP_RART_NUMERIC = UInt32( 138 )
global const MATOP_SET_BLOCK_SIZES = UInt32( 139 )
global const MATOP_AYPX = UInt32( 140 )
global const MATOP_RESIDUAL = UInt32( 141 )
global const MATOP_FDCOLORING_SETUP = UInt32( 142 )
global const MATOP_MPICONCATENATESEQ = UInt32( 144 )

global const InsertMode = Cint
global const INSERT_VALUES = Cint( 1 )
global const ADD_VALUES = Cint( 2 )

global const KSPType = String
global const KSPRICHARDSON = String( "richardson" )
global const KSPCHEBYSHEV = String( "chebyshev" )
global const KSPCG = String( "cg" )
global const KSPGROPPCG = String( "groppcg" )
global const KSPPIPECG = String( "pipecg" )
global const KSPCGNE = String( "cgne" )
global const KSPNASH = String( "nash" )
global const KSPSTCG = String( "stcg" )
global const KSPGLTR = String( "gltr" )
global const KSPFCG = String( "fcg" )
global const KSPGMRES = String( "gmres" )
global const KSPFGMRES = String( "fgmres" )
global const KSPLGMRES = String( "lgmres" )
global const KSPDGMRES = String( "dgmres" )
global const KSPPGMRES = String( "pgmres" )
global const KSPTCQMR = String( "tcqmr" )
global const KSPBCGS = String( "bcgs" )
global const KSPIBCGS = String( "ibcgs" )
global const KSPFBCGS = String( "fbcgs" )
global const KSPFBCGSR = String( "fbcgsr" )
global const KSPBCGSL = String( "bcgsl" )
global const KSPCGS = String( "cgs" )
global const KSPTFQMR = String( "tfqmr" )
global const KSPCR = String( "cr" )
global const KSPPIPECR = String( "pipecr" )
global const KSPLSQR = String( "lsqr" )
global const KSPPREONLY = String( "preonly" )
global const KSPQCG = String( "qcg" )
global const KSPBICG = String( "bicg" )
global const KSPMINRES = String( "minres" )
global const KSPSYMMLQ = String( "symmlq" )
global const KSPLCD = String( "lcd" )
global const KSPPYTHON = String( "python" )
global const KSPGCR = String( "gcr" )

global const PCJacobiType = Int32
global const PC_JACOBI_DIAGONAL = Int32( 0 )
global const PC_JACOBI_ROWMAX = Int32( 1 )
global const PC_JACOBI_ROWSUM = Int32( 2 )

global const PCType = String
global const PCNONE = String( "none" )
global const PCJACOBI = String( "jacobi" )
global const PCSOR = String( "sor" )
global const PCLU = String( "lu" )
global const PCSHELL = String( "shell" )
global const PCBJACOBI = String( "bjacobi" )
global const PCMG = String( "mg" )
global const PCEISENSTAT = String( "eisenstat" )
global const PCILU = String( "ilu" )
global const PCICC = String( "icc" )
global const PCASM = String( "asm" )
global const PCGASM = String( "gasm" )
global const PCKSP = String( "ksp" )
global const PCCOMPOSITE = String( "composite" )
global const PCREDUNDANT = String( "redundant" )
global const PCSPAI = String( "spai" )
global const PCNN = String( "nn" )
global const PCCHOLESKY = String( "cholesky" )
global const PCPBJACOBI = String( "pbjacobi" )
global const PCMAT = String( "mat" )
global const PCHYPRE = String( "hypre" )
global const PCPARMS = String( "parms" )
global const PCFIELDSPLIT = String( "fieldsplit" )
global const PCTFS = String( "tfs" )
global const PCML = String( "ml" )
global const PCGALERKIN = String( "galerkin" )
global const PCEXOTIC = String( "exotic" )
global const PCCP = String( "cp" )
global const PCBFBT = String( "bfbt" )
global const PCLSC = String( "lsc" )
global const PCPYTHON = String( "python" )
global const PCPFMG = String( "pfmg" )
global const PCSYSPFMG = String( "syspfmg" )
global const PCREDISTRIBUTE = String( "redistribute" )
global const PCSVD = String( "svd" )
global const PCGAMG = String( "gamg" )
global const PCSACUSP = String( "sacusp" )
global const PCSACUSPPOLY = String( "sacusppoly" )
global const PCBICGSTABCUSP = String( "bicgstabcusp" )
global const PCAINVCUSP = String( "ainvcusp" )
global const PCBDDC = String( "bddc" )
global const PCKACZMARZ = String( "kaczmarz" )

global const MatAssemblyType = Cint
global const MAT_FLUSH_ASSEMBLY = Cint( 1 )
global const MAT_FINAL_ASSEMBLY = Cint( 0 )

global const VecType = String
global const VECSEQ = String( "seq" )
global const VECMPI = String( "mpi" )
global const VECSTANDARD = String( "standard" )
global const VECSHARED = String( "shared" )
global const VECSEQCUSP = String( "seqcusp" )
global const VECMPICUSP = String( "mpicusp" )
global const VECCUSP = String( "cusp" )
global const VECSEQVIENNACL = String( "seqviennacl" )
global const VECMPIVIENNACL = String( "mpiviennacl" )
global const VECVIENNACL = String( "viennacl" )
global const VECSEQCUDA = String( "seqcuda" )
global const VECMPICUDA = String( "mpicuda" )
global const VECCUDA = String( "cuda" )
global const VECNEST = String( "nest" )

global const PCSide = Cint
global const PC_SIDE_DEFAULT = Cint( -1 )
global const PC_LEFT = Cint( 0 )
global const PC_RIGHT = Cint( 1 )
global const PC_SYMMETRIC = Cint( 2 )

global const MatType = String
global const MATSAME = String( "same" )
global const MATMAIJ = String( "maij" )
global const MATSEQMAIJ = String( "seqmaij" )
global const MATMPIMAIJ = String( "mpimaij" )
global const MATIS = String( "is" )
global const MATAIJ = String( "aij" )
global const MATSEQAIJ = String( "seqaij" )
global const MATMPIAIJ = String( "mpiaij" )
global const MATAIJCRL = String( "aijcrl" )
global const MATSEQAIJCRL = String( "seqaijcrl" )
global const MATMPIAIJCRL = String( "mpiaijcrl" )
global const MATAIJCUSP = String( "aijcusp" )
global const MATSEQAIJCUSP = String( "seqaijcusp" )
global const MATMPIAIJCUSP = String( "mpiaijcusp" )
global const MATAIJCUSPARSE = String( "aijcusparse" )
global const MATSEQAIJCUSPARSE = String( "seqaijcusparse" )
global const MATMPIAIJCUSPARSE = String( "mpiaijcusparse" )
global const MATAIJVIENNACL = String( "aijviennacl" )
global const MATSEQAIJVIENNACL = String( "seqaijviennacl" )
global const MATMPIAIJVIENNACL = String( "mpiaijviennacl" )
global const MATAIJPERM = String( "aijperm" )
global const MATSEQAIJPERM = String( "seqaijperm" )
global const MATMPIAIJPERM = String( "mpiaijperm" )
global const MATSHELL = String( "shell" )
global const MATDENSE = String( "dense" )
global const MATSEQDENSE = String( "seqdense" )
global const MATMPIDENSE = String( "mpidense" )
global const MATELEMENTAL = String( "elemental" )
global const MATBAIJ = String( "baij" )
global const MATSEQBAIJ = String( "seqbaij" )
global const MATMPIBAIJ = String( "mpibaij" )
global const MATMPIADJ = String( "mpiadj" )
global const MATSBAIJ = String( "sbaij" )
global const MATSEQSBAIJ = String( "seqsbaij" )
global const MATMPISBAIJ = String( "mpisbaij" )
global const MATSEQBSTRM = String( "seqbstrm" )
global const MATMPIBSTRM = String( "mpibstrm" )
global const MATBSTRM = String( "bstrm" )
global const MATSEQSBSTRM = String( "seqsbstrm" )
global const MATMPISBSTRM = String( "mpisbstrm" )
global const MATSBSTRM = String( "sbstrm" )
global const MATDAAD = String( "daad" )
global const MATMFFD = String( "mffd" )
global const MATNORMAL = String( "normal" )
global const MATLRC = String( "lrc" )
global const MATSCATTER = String( "scatter" )
global const MATBLOCKMAT = String( "blockmat" )
global const MATCOMPOSITE = String( "composite" )
global const MATFFT = String( "fft" )
global const MATFFTW = String( "fftw" )
global const MATSEQCUFFT = String( "seqcufft" )
global const MATTRANSPOSEMAT = String( "transpose" )
global const MATSCHURCOMPLEMENT = String( "schurcomplement" )
global const MATPYTHON = String( "python" )
global const MATHYPRESTRUCT = String( "hyprestruct" )
global const MATHYPRESSTRUCT = String( "hypresstruct" )
global const MATSUBMATRIX = String( "submatrix" )
global const MATLOCALREF = String( "localref" )
global const MATNEST = String( "nest" )

global const NormType = Int32
global const NORM_1 = Int32( 0 )
global const NORM_2 = Int32( 1 )
global const NORM_FROBENIUS = Int32( 2 )
global const NORM_INFINITY = Int32( 3 )
global const NORM_MAX = Int32( 3 )

global const PetscMatStructure = Int32
global const DIFFERENT_NONZERO_PATTERN = Int32( 0 )
global const SUBSET_NONZERO_PATTERN = Int32( 1 )
global const SAME_NONZERO_PATTERN = Int32( 2 )

global const MatReuse = Int32
global const MAT_INITIAL_MATRIX = Int32( 0 )
global const MAT_REUSE_MATRIX = Int32( 1 )
global const MAT_IGNORE_MATRIX = Int32( 2 )
global const MAT_INPLACE_MATRIX = Int32( 3 )


text:
struct controls
    ref_changes::Dict{Any, Any}
    load_changes::Dict{Any, Any}
    line_trips::Dict{Any, Any}
    source_trip::Dict{Any, Any}
    line_sc::Dict{Any, Any}
end


text:
module RulesBernoulliOutTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_rules

@testset "rules:Bernoulli:out" begin

    @testset "Belief Propagation: (m_p::PointMass)" begin

        @test_rules [ with_float_conversions = true ] Bernoulli(:out, Marginalisation) [
            (input = (m_p = PointMass(1.0),), output = Bernoulli(1.0)),
            (input = (m_p = PointMass(0.2),), output = Bernoulli(0.2)),
        ]

    end

    @testset "Variational Message Passing: (q_p::PointMass)" begin

        @test_rules [ with_float_conversions = true ] Bernoulli(:out, Marginalisation) [
            (input = (q_p = PointMass(1.0),), output = Bernoulli(1.0)),
            (input = (q_p = PointMass(0.3),), output = Bernoulli(0.3)),
        ]

    end

    @testset "Variational Message Passing: (q_p::Beta)" begin

        @test_rules [ with_float_conversions = true ] Bernoulli(:out, Marginalisation) [
            (input = (q_p = Beta(1.0, 1.0),), output = Bernoulli(0.5)),
            (input = (q_p = Beta(0.2, 0.2),), output = Bernoulli(0.5)),
        ]

    end
end
end


text:
iscancread(ncname::AbstractString,fol::AbstractString="") = Dataset(joinpath(fol,ncname))

function iscarawfolder(iroot::AbstractDict;irun::Integer)

    return joinpath(iroot["raw"],"run$(@sprintf("%04d",irun))")

end

function iscarawname(iroot::AbstractDict;irun::Integer)

    return joinpath(iroot["raw"],"run$(@sprintf("%04d",irun))",iroot["fname"])

end

function iscarawread(ipar::AbstractDict,iroot::AbstractDict;irun::Integer)

    inc = iscarawname(iroot;irun=irun); ids = iscancread(inc)
    return ids,ids[ipar["ID"]]

end

function iscacalcname(ipar::AbstractDict,iroot::AbstractDict;irun::Integer)

    return joinpath(iroot["raw"],"run$(@sprintf("%04d",irun))","atmos_$(ipar["ID"]).nc")

end

function iscacalcread(ipar::AbstractDict,iroot::AbstractDict;irun::Integer)

    inc = iscacalcname(ipar,iroot;irun=irun); ids = iscancread(inc)
    return ids,ids[ipar["ID"]]

end

function iscaanafolder(ipar::AbstractDict,iroot::AbstractDict)

    if ipar["level"] == "sfc";

        fol = joinpath(iroot["ana"],ipar["ID"]);
        if !isdir(fol)
            @info "$(Dates.now()) - The folder for analyzed $(uppercase(ipar["ID"])) data does not exist.  Creating now ..."
            mkpath(fol);
        end

    else

        phPa = "$(ipar["ID"])-lvl$(@sprintf("%02d",ipar["level"]))"
        fol = joinpath(iroot["ana"],ipar["ID"],phPa);
        if !isdir(fol)
            @info "$(Dates.now()) - The folder for analyzed $(uppercase(ipar["ID"])) data at σ-level $(ipar["level"]) does not exist.  Creating now ..."
            mkpath(fol);
        end

    end

    return fol
end

function iscaananame(ipar::AbstractDict;irun::Integer)

    if ipar["level"] != "sfc"
          fname = "$(ipar["ID"])-lvl$(@sprintf("%02d",ipar["level"]))";
    else; fname = "$(ipar["ID"])-sfc";
    end

    return "$(fname)-run$(@sprintf("%04d",irun)).nc"

end

function iscaanaread(
    ID::AbstractString,
    ipar::AbstractDict,iroot::AbstractDict;
    irun::Integer
)

    ibase = iscaanafol(ipar,iroot)
    inc = iscaananame(ipar;irun=irun);
    ids = iscancread(inc,ibase)
    return ids,ids[ID]

end

function iscapre2lvl(pressure::Real,imod::AbstractDict)

    return argmin(abs.(imod["levels"] .- pressure))

end

function putinfo(imod::Dict,ipar::Dict,itime::Dict,iroot::Dict)

    rfol = pwd();
    cd(iscaanafolder(ipar,iroot)); @save "info_par.jld2" imod ipar;
    cd(iroot["ana"]); @save "info_time.jld2" itime;
    cd(rfol);

end


text:
#
# Autogenerated by Thrift Compiler (0.11.0)
#
# DO NOT EDIT UNLESS YOU ARE SURE THAT YOU KNOW WHAT YOU ARE DOING

struct _enum_TResultSetLayout
  GroupByPerfectHash::Int32
  GroupByBaselineHash::Int32
  Projection::Int32
  NonGroupedAggregate::Int32
end
const TResultSetLayout = _enum_TResultSetLayout(Int32(0), Int32(1), Int32(2), Int32(3))

struct _enum_TCountDistinctImplType
  Invalid::Int32
  Bitmap::Int32
  StdSet::Int32
end
const TCountDistinctImplType = _enum_TCountDistinctImplType(Int32(0), Int32(1), Int32(2))

struct _enum_TAggKind
  AVG::Int32
  MIN::Int32
  MAX::Int32
  SUM::Int32
  COUNT::Int32
  APPROX_COUNT_DISTINCT::Int32
  SAMPLE::Int32
end
const TAggKind = _enum_TAggKind(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6))

mutable struct TSlotSize <: Thrift.TMsg
  padded::Int16
  logical::Int16
  TSlotSize() = (o=new(); fillunset(o); o)
end # mutable struct TSlotSize

mutable struct TColSlotContext <: Thrift.TMsg
  slot_sizes::Vector{TSlotSize}
  col_to_slot_map::Vector{Vector{Int32}}
  TColSlotContext() = (o=new(); fillunset(o); o)
end # mutable struct TColSlotContext

mutable struct TCountDistinctDescriptor <: Thrift.TMsg
  impl_type::Int32
  min_val::Int64
  bitmap_sz_bits::Int64
  approximate::Bool
  device_type::Int32
  sub_bitmap_count::Int64
  TCountDistinctDescriptor() = (o=new(); fillunset(o); o)
end # mutable struct TCountDistinctDescriptor

const TCountDistinctDescriptors = Vector{TCountDistinctDescriptor}  #moved this to end to be after TCountDistinctDescriptor

mutable struct TResultSetBufferDescriptor <: Thrift.TMsg
  layout::Int32
  keyless::Bool
  entry_count::Int32
  idx_target_as_key::Int32
  min_val::Int64
  max_val::Int64
  bucket::Int64
  group_col_widths::Vector{Int16}
  key_bytewidth::Int16
  col_slot_context::TColSlotContext
  target_groupby_indices::Vector{Int32}
  count_distinct_descriptors::TCountDistinctDescriptors
  force_4byte_float::Bool
  TResultSetBufferDescriptor() = (o=new(); fillunset(o); o)
end # mutable struct TResultSetBufferDescriptor

mutable struct TTargetInfo <: Thrift.TMsg
  is_agg::Bool
  kind::Int32
  _type::TTypeInfo
  arg_type::TTypeInfo
  skip_nulls::Bool
  is_distinct::Bool
  TTargetInfo() = (o=new(); fillunset(o); o)
end # mutable struct TTargetInfo

mutable struct TCountDistinctSetStorage <: Thrift.TMsg
  bitmap::Vector{UInt8}
  sparse_set::Set{Int64}
  TCountDistinctSetStorage() = (o=new(); fillunset(o); o)
end # mutable struct TCountDistinctSetStorage
meta(t::Type{TCountDistinctSetStorage}) = meta(t, Symbol[:bitmap,:sparse_set], Int[], Dict{Symbol,Any}())

mutable struct TCountDistinctSet <: Thrift.TMsg
  _type::Int32
  storage::TCountDistinctSetStorage
  remote_ptr::Int64
  TCountDistinctSet() = (o=new(); fillunset(o); o)
end # mutable struct TCountDistinctSet

mutable struct TSerializedRows <: Thrift.TMsg
  buffers::Vector{Vector{UInt8}}
  buffer_lengths::Vector{Int64}
  buffers_total_size::Int64
  total_compression_time_ms::Int32
  descriptor::TResultSetBufferDescriptor
  targets::Vector{TTargetInfo}
  target_init_vals::Vector{Int64}
  varlen_buffer::Vector{Vector{UInt8}}
  count_distinct_sets::Vector{TCountDistinctSet}
  explanation::String
  TSerializedRows() = (o=new(); fillunset(o); o)
end # mutable struct TSerializedRows


text:
# lists table-specific iterators
iters(session::AccumuloSession, tablename::String) = listIterators(client(session), handle(session), tablename)

# get a single iterator with its settings
iter(session::AccumuloSession, tablename::String, itername::String, scope::Integer) = getIteratorSetting(client(session), handle(session), tablename, itername, Int32(scope))

# create an iterator type
function iter(name::String, class::String, priority::Integer, properties::Dict=Dict())
    iter_properties = Dict{String,String}(string(n)=>string(v) for (n,v) in properties)

    thriftbuild(IteratorSetting, Dict(:priority => Int32(priority),
        :name => name,
        :iteratorClass => class,
        :properties => iter_properties))
end

# attach an iterator to a table
function add_iter(session::AccumuloSession, tablename::String, iter_setting::IteratorSetting, scopes::SET; check::Bool=true)
    iter_scopes = Set([Int32(scope) for scope in scopes])

    if check
        checkIteratorConflicts(client(session), handle(session), tablename, iter_setting, iter_scopes)
    end
    attachIterator(client(session), handle(session), tablename, iter_setting, iter_scopes)
end

# detach the iterator from the table
remove_iter(session::AccumuloSession, tablename::String, iter_setting::IteratorSetting, scopes::SET) = remove_iter(session, tablename, iter_setting.name, scopes)
function remove_iter(session::AccumuloSession, tablename::String, name::String, scopes::SET)
    iter_scopes = Set([Int32(scope) for scope in scopes])
    removeIterator(client(session), handle(session), tablename, iter_scopes)
end

# check if there is any conflicts with the iterator name and priority
function check_iter(session::AccumuloSession, tablename::String, iter_setting::IteratorSetting, scopes::SET)
    iter_scopes = Set([Int32(scope) for scope in scopes])
    checkIteratorConflicts(client(session), handle(session), tablename, iter_setting, iter_scopes)
end


text:
using Jchemo

n = 100 ; p = 10
X = rand(n, p)
y = rand([3, 4, 10], n)
#y = rand(["a", "b", "c"], n) ;
Xtrain = X[1:70, :] ; ytrain = y[1:70]
Xtest = X[71:n, :] ; ytest = y[71:n]

tab(ytrain)
tab(ytest)

############ PLSR-DA

nlv = 0:3
pars = mpars(nlv = nlv)
gridscore(Xtrain, ytrain, Xtest, ytest;
    score = err, fun = plsrda, pars = pars)

# Faster
gridscorelv(Xtrain, ytrain, Xtest, ytest;
    score = err, fun = plsrda, nlv = nlv)

############ PLS-LDA/QDA

nlv = 1:3
pars = mpars(nlv = nlv)
gridscore(Xtrain, ytrain, Xtest, ytest;
    score = err, fun = plslda, pars = pars)

# Faster
gridscorelv(Xtrain, ytrain, Xtest, ytest;
    score = err, fun = plslda, nlv = nlv)

gridscorelv(Xtrain, ytrain, Xtest, ytest;
    score = err, fun = plsqda, nlv = nlv)

############ RR-DA

lb = [.001; .01; .1]
pars = mpars(lb = lb)
gridscore(Xtrain, ytrain, Xtest, ytest;
    score = err, fun = rrda, pars = pars)

# Faster
gridscorelb(Xtrain, ytrain, Xtest, ytest;
    score = err, fun = rrda, lb = lb)

############ KRR-DA

lb = [.001; .01; .1]
pars = mpars(lb = lb, gamma = [1; 10])
#pars = mpars(lb = lb, kern = ["kpol";], degree = 1:2)
gridscore(Xtrain, ytrain, Xtest, ytest;
    score = err, fun = krrda, pars = pars)

# Faster
pars = mpars(gamma = [1; 10])
gridscorelb(Xtrain, ytrain, Xtest, ytest;
    score = err, fun = krrda, lb = lb, pars = pars)

############ kNN-DA

nlvdis = 5 ; metric = ["mahal";]
h = [1; 3] ; k = [20; 10]
pars = mpars(nlvdis = nlvdis, metric = metric, h = h, k = k)
gridscore(Xtrain, ytrain, Xtest, ytest;
    score = msep, fun = knnda, pars = pars, verbose = true)

############ kNN-LWPLSR-DA

nlvdis = 5 ; metric = ["mahal";]
h = [1; 3] ; k = [20; 10]
nlv = 1:2
pars = mpars(nlv = nlv, nlvdis = nlvdis, metric = metric, h = h, k = k)
gridscore(Xtrain, ytrain, Xtest, ytest;
    score = msep, fun = lwplsrda, pars = pars, verbose = true)

# Faster
pars = mpars(nlvdis = nlvdis, metric = metric, h = h, k = k)
gridscorelv(Xtrain, ytrain, Xtest, ytest;
    score = msep, fun = lwplsrda, nlv = nlv, pars = pars, verbose = false)

############ kNN-LWPLS-LDA

nlvdis = 5 ; metric = ["mahal";]
h = [1; 3] ; k = [20; 10]
nlv = 1:2
pars = mpars(nlv = nlv, nlvdis = nlvdis, metric = metric, h = h, k = k)
gridscore(Xtrain, ytrain, Xtest, ytest;
    score = msep, fun = lwplslda, pars = pars, verbose = true)

# Faster
pars = mpars(nlvdis = nlvdis, metric = metric, h = h, k = k)
gridscorelv(Xtrain, ytrain, Xtest, ytest;
    score = msep, fun = lwplslda, nlv = nlv, pars = pars, verbose = false)

########### PLSR-DA-AGG

# Here there is no sense to use gridscorelv

pars = mpars(nlv = ["1:2"; "1:3"])
gridscore(Xtrain, ytrain, Xtest, ytest;
    score = msep, fun = plsrda_agg, pars = pars)

############ kNN-LWPLSR-DA-AGG

# Here there is no sense to use gridscorelv

nlvdis = 5 ; metric = ["mahal";] ;
h = [1.; 3.] ; k = [20; 10]
nlv = ["1:2"; "2:5"] ;
pars = mpars(nlv = nlv, nlvdis = nlvdis, metric = metric, h = h, k = k) ;
gridscore(Xtrain, ytrain, Xtest, ytest;
    score = msep, fun = lwplsrda_agg, pars = pars, verbose = true)


text:
#
# This file contains various distributions in unconstrained, centered, and scaled form.
# Sometimes ad hoc, via autodiff.
# The purpose is for finding error-minimizing (Optimal) quadrature rules.
#

## Normal

## Gamma - low var

## Gamma - medium var

## Gamma - high var

## Inv-Gamma - low var

## Inv-Gamma - medium var

## Inv-Gamma - high var

## T - small nu

## T - medium nu

## T - large nu


text:
"""
    PrettyLogging

A module to help simplify "verbose" output from complicated functions/scripts.

Goal for usage:
```julia
function complicatedfun(a, b; verbose=0, kwargs...)
    # do some preprocessing
    @logmsg "preprocessing complete"
    # now call a subroutine
    @logmsg foo(b)

end
```
"""
module PrettyLogging

const INDENT_WIDTH = 2
const INDENT = "  "
const START = "<<< "
const END = ">>> "
const IMPORTANT = "* "
const ERROR = "! "
const DATE_FMT = "YYYY mm dd HH:MM:SS"

"""
    Logger

Object that handles writing log messages
"""
mutable struct Logger
    io::IO
    disp::Int # display threshold; all messages get flushed to io but not necessarily printed to stdout
    hold
    indent::Int
end
Logger(filename="$(@__FILE__).$(myid()).log"; disp=0) = Logger(open(filename, "a"), disp, "", 1)

"""
    indent(L::Logger)

Returns a new `Logger` with the same IO stream and display settings and indentation
increased by 1. Makes it easy to have output from nested function calls indented, if
desired. Clears the `hold` field.

"""
indent(L::Logger) = Logger(L.io, L.disp, "", L.indent+1)

"""
    logmsg(L::logger, m, v=0; timestamp=true)

Write message `m` to the log file (i.e. `L.io`). If `L` is holding any messages, flush that
message to the log file as well and clear the hold. If `v <= L.disp`, print the output to
`STDOUT` as well.

"""
function logmsg(L::Logger, m, v=0; note=INDENT)
    flushhold(L)
    msg = "$(ts()) $note$(spaces(L.indent-1))$m"
    println(L.io, msg)
    if v <= L.disp
        println(msg)
    end
end
logerr(L::Logger, m, v=0) = logmsg(L, m, v; note=ERROR)
logspecial(L::Logger, m, v=0) = logmsg(L, m, v; note=IMPORTANT)

function loghold(L::Logger, m; extra_space=0)
    if length(L.hold) == 0
        L.hold = ts() * " " * spaces(L.indent + extra_space)
    end
    L.hold = L.hold * m
end

function flushhold(L::Logger)
    if length(L.hold) > 0
        println(L.io, L.hold)
        L.hold = ""
    end
end

"""
    ts()

Returns the current time as a formatted string for logging purposes.
"""
ts(fmt=DATE_FMT) = Dates.format(now(), fmt)

spaces(n=1) = INDENT ^ n

export Logger, indent,
       logmsg, logerr, logspecial, loghold,
       ts

end # module


text:
# ==============================================================================
# The gravity machine (Man of Steel) -> to terraform the world

println("""\nAlgorithme "Gravity machine" --------------------------------\n""")

const verbose = true
const graphic = true

println("-) Active les packages requis\n")
using JuMP, GLPK, PyPlot, Printf, Random
verbose ? println("  Fait \n") : nothing

include("GMparsers.jl")        # parsers of instances and non-dominated points
include("GMjumpModels.jl")     # JuMP models for computing optima
include("GMmopPrimitives.jl")  # usuals algorithms in multiobjective optimization

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

# types ------------------------------------------------------------------------

# type corresponding to a solution
mutable struct tSolution{T}
    x :: Vector{T}                # vector variables x (1..n)
    y :: Vector{T}                # vector outcomes  y (1..p)
end

# type corresponding to a point generator
mutable struct tGenerateur
    sRel :: tSolution{Float64}    # initial relaxed solution
    sInt :: tSolution{Int64}      # integer solution
    sPrj :: tSolution{Float64}    # projected solution
    sFea :: Bool                  # indicate if sInt is feasible or not
end

# type of a point (x,y) in the objective space
mutable struct tPoint
    x :: Float64
    y :: Float64
end

# Global variables -------------------------------------------------------------

# listes de points pour les affichages graphiques
xLf1  = (Float64)[]; yLf1  = (Float64)[] # liste des points (x,y) relaches
xLf2  = (Float64)[]; yLf2  = (Float64)[] # liste des points (x,y) relaches
xL    = (Float64)[]; yL    = (Float64)[] # liste des points (x,y) relaches
XInt  = (Int64)[];   YInt  = (Int64)[]   # liste des points (x,y) entiers
XProj = (Float64)[]; YProj = (Float64)[] # liste des points (x,y) projetes
XFeas = (Int64)[];   YFeas = (Int64)[]   # liste des points (x,y) admissibles
XPert = (Int64)[];   YPert = (Int64)[]   # liste des points (x,y) perturbes

# ==============================================================================
# Initialisation structure donnees contenant tous les generateurs

function allocateDatastructure(nbgen::Int64, nbvar::Int64, nbobj::Int64, )

    verbose ? println("\n  → Allocation memoire pour ",nbgen," generateurs\n") : nothing

    vg = Vector{tGenerateur}(undef, nbgen)
    for k = 1:nbgen
        vg[k] = tGenerateur( tSolution{Float64}(zeros(Float64,nbvar),zeros(Float64,nbobj)),
                              tSolution{Int64}(zeros(Int64,nbvar),zeros(Int64,nbobj)),
                              tSolution{Float64}(zeros(Float64,nbvar),zeros(Float64,nbobj)),
                              false
                            )
    end
    return vg
end

# ==============================================================================
# Ajout d'une solution relachee initiale a un generateur

function ajouterX0!(vg::Vector{tGenerateur}, k::Int64, s::tSolution{Float64})

    vg[k].sRel = deepcopy(s) # met en place le generateur \bar{x}^k
    vg[k].sPrj = deepcopy(s) # le generateur est la premiere projection \bar{x}^{k,0}
    return nothing
end

# ==============================================================================
# Ajout d'une solution entiere (arrondie ou perturbee) a un generateur

function ajouterXtilde!(vg::Vector{tGenerateur}, k::Int64, x::Vector{Int64}, y::Vector{Int64})

    vg[k].sInt.x = copy(x)
    vg[k].sInt.y = copy(y)
    return nothing
end

# ==============================================================================
# Ajout d'une solution fractionnaire (projetee) a un generateur

function ajouterXbar!(vg::Vector{tGenerateur}, k::Int64, x::Vector{Float64}, y::Vector{Float64})

    vg[k].sPrj.x = copy(x)
    vg[k].sPrj.y = copy(y)
    return nothing
end

# ==============================================================================
# Elabore 2 ensembles d'indices selon que xTilde[i] vaut 0 ou 1

function split01(xTilde)

   indices0 = (Int64)[]
   indices1 = (Int64)[]

   for i=1:length(xTilde)
       if xTilde[i] == 0
           push!(indices0,i)
       else
           push!(indices1,i)
       end
    end

   return indices0, indices1
end

# ==============================================================================
# Projete xTilde sur le polyedre X

function Δ2SPA(A::Array{Int,2}, xTilde::Array{Int,1})

    nbctr = size(A,1)
    nbvar = size(A,2)
    idxTilde0, idxTilde1 = split01(xTilde)

    proj = Model(GLPK.Optimizer)
    @variable(proj, 0.0 <= x[1:length(xTilde)] <= 1.0 )
    @objective(proj, Min, sum(x[i] for i in idxTilde0) + sum((1-x[i]) for i in idxTilde1) )
    @constraint(proj, [i=1:nbctr],(sum((x[j]*A[i,j]) for j in 1:nbvar)) == 1)
    optimize!(proj)
    return objective_value(proj), value.(x)
end

# ==============================================================================
# Projete xTilde sur le polyedre X

function Δ2SPAbis(A::Array{Int,2}, xTilde::Array{Int,1}, c1, c2, k, λ1, λ2)

    nbctr = size(A,1)
    nbvar = size(A,2)
    idxTilde0, idxTilde1 = split01(xTilde)

#    cλ = 0.5.*c1 + 0.5.*c2
    proj = Model(GLPK.Optimizer)
    @variable(proj, 0.0 <= x[1:length(xTilde)] <= 1.0 )
    @objective(proj, Min, sum(λ1[k]*x[i] for i in idxTilde0) + sum(λ2[k]*(1-x[i]) for i in idxTilde1) )
#    @objective(proj, Min, sum(cλ[i]*x[i] for i in idxTilde0) + sum(cλ[i]*(1-x[i]) for i in idxTilde1) )
    @constraint(proj, [i=1:nbctr],(sum((x[j]*A[i,j]) for j in 1:nbvar)) == 1)
    optimize!(proj)
    return objective_value(proj), value.(x)
end

# ==============================================================================
# test si une solution est admissible en verifiant si sa relaxation lineaire
# conduit a une solution entiere

function estAdmissible(x)

    admissible = true
    i=1
    while admissible && i<=length(x)
        if round(x[i], digits=3)!=0.0 && round(x[i], digits=3)!=1.0
            admissible = false
        end
        i+=1
    end
    return admissible
end

# ==============================================================================
# calcule la performance z d'une solution x sur les 2 objectifs

function evaluerSolution(x, c1, c2)

    z1 = 0.0; z2 = 0.0
    for i in 1:length(x)
        z1 += x[i] * c1[i]
        z2 += x[i] * c2[i]
    end
    return round(z1, digits=2), round(z2, digits=2)
end

# ==============================================================================
# nettoyage des valeurs des variables d'une solution x relachee sur [0,1]

function nettoyageSolution!(x)

    nbvar=length(x)
    for i in 1:nbvar
        if     round(x[i], digits=3) == 0.0
                   x[i] = 0.0
        elseif round(x[i], digits=3) == 1.0
                   x[i] = 1.0
        else
                   x[i] = round(x[i], digits=3)
        end
    end
end

# ==============================================================================
# Calcul des generateurs avec 2 ϵ-contraintes alternees jusqu'a leur rencontre

function calculGenerateurs(A, c1, c2, tailleSampling, minf1RL, maxf2RL, maxf1RL, minf2RL)

    # acces aux var globales : xLf1, yLf1, xLf2, yLf2

    nbctr = size(A,1)
    nbvar = size(A,2)

    L1 = (tSolution{Float64})[]
    L2 = (tSolution{Float64})[]

    # Premiere generation de points (avec z1 la fonction a minimiser) ----------
    pasSample2 = (maxf2RL - minf2RL) / (tailleSampling-1) # pas de l'echantillonage sur z2
    j1 = 1

    # Seconde generation de points (avec z2 la fonction a minimiser) -----------
    pasSample1 = (maxf1RL-minf1RL) / (tailleSampling-1) # pas de l'echantillonage sur z1
    j2 = 1

    alternance = 1
    maxf2RLlimite = maxf2RL
    minf2RLlimite = minf2RL

    while maxf2RLlimite > minf2RLlimite
        alternance += 1
        if alternance % 2 == 0

            # minimise f1 avec ϵ-contrainte sur f2 -----------------------------
            verbose ? @printf("  z1 %2d : ϵ = %8.2f  ", j1, maxf2RL - (j1-1) * pasSample2) : nothing # echantillonage sur z2

            # calcul d'une solution epsilon-contrainte
            f1RL, xf1RL = computeLinearRelax2SPA(nbvar, nbctr, A, c1, c2, maxf2RL - (j1-1) * pasSample2, 1)

            # reconditionne les valeurs 0 et 1 et arrondi les autres valeurs
            nettoyageSolution!(xf1RL)
            verbose ? @printf("fRL = %8.2f  ",round(f1RL, digits=2)) : nothing

            # recalcule la solution au regard des 2 objectifs
            z1f1RLcourant, z2f1RLcourant = evaluerSolution(xf1RL, c1, c2)
            verbose ? @printf("[ %8.2f , %8.2f ] \n", z1f1RLcourant, z2f1RLcourant) : nothing

            # maj la valeur limite sur l'objectif 2 pour la solution courante
            maxf2RLlimite = z2f1RLcourant
            if maxf2RLlimite > minf2RLlimite
                push!(L1, (tSolution{Float64})(xf1RL,[z1f1RLcourant, z2f1RLcourant]))
                push!(xL,z1f1RLcourant);push!(yL,z2f1RLcourant)
                push!(xLf1,z1f1RLcourant);push!(yLf1,z2f1RLcourant)
            end
            j1 = j1+1

        else

            # minimise f2 avec ϵ-contrainte sur f1 -----------------------------
            verbose ? @printf("  z2 %2d : ϵ = %8.2f  ", j2, maxf1RL - (j2-1) * pasSample1) : nothing # echantillonage sur z1

            # calcul d'une solution epsilon-contrainte
            f2RL, xf2RL = computeLinearRelax2SPA(nbvar, nbctr, A, c1, c2, maxf1RL - (j2-1) * pasSample1, 2)

            # reconditionne les valeurs 0 et 1 et arrondi les autres valeurs
            nettoyageSolution!(xf2RL)
            verbose ? @printf("fRL = %8.2f  ",round(f2RL, digits=2)) : nothing

            # recalcule la solution au regard des 2 objectifs
            z1f2RLcourant, z2f2RLcourant = evaluerSolution(xf2RL, c1, c2)
            verbose ? @printf("[ %8.2f , %8.2f ] \n", z1f2RLcourant, z2f2RLcourant) : nothing

            # maj la valeur limite sur l'objectif 2 pour la solution courante
            minf2RLlimite = z2f2RLcourant
            if maxf2RLlimite > minf2RLlimite
                push!(L2, (tSolution{Float64})(xf2RL,[z1f2RLcourant, z2f2RLcourant]))
                push!(xL,z1f2RLcourant);push!(yL,z2f2RLcourant)
                push!(xLf2,z1f2RLcourant);push!(yLf2,z2f2RLcourant)
            end
            j2 = j2+1
        end
    end
    return  length(L1)+length(L2), vcat(L1,reverse(L2))
end

# ==============================================================================
# predicat : verifie si une solution entiere est realisable
function isFeasible(vg,k)
    #verbose && vg[k].sFea == true ? println("   feasible") : nothing
    return (vg[k].sFea == true)
end

# ==============================================================================
# predicat : verifie si le nombre d'essai maximum a ete tente
function isFinished(trial, maxTrial)
#    verbose && trial > maxTrial ? println("   maxTrial") : nothing
    return (trial > maxTrial)
end

# ==============================================================================
# predicat : verifie si le budget de calcul maximum a ete consomme
function isTimeout(temps, maxTime)
#    verbose && time()- temps > maxTime ? println("   maxTime") : nothing
    return (time()- temps > maxTime)
end

# ==============================================================================
# elabore pC le pointeur du cone ouvert vers L

function elaborePointConeOuvertversL(vg, k, pB, pA)

    # recupere les coordonnees du point projete
    pC=tPoint(vg[k].sPrj.y[1], vg[k].sPrj.y[2])

    # etablit le point nadir pN au depart des points pA et pB adjacents au generateur k
    pN = tPoint( pA.x , pB.y )

#    print("Coordonnees du cone 2 : ")
#    @show pC, pN

    # retient pN si pC domine pN (afin d'ouvrir le cone)
    if pC.x < pN.x  &&  pC.y < pN.y
        # remplace pC par pN
        pC=tPoint( pA.x , pB.y )
    end

    return pC
end

# ==============================================================================
#= Retourne un booléen indiquant si un point se trouve dans un secteur défini dans
  le sens de rotation trigonométrique (repère X de gauche à droite, Y du haut vers
  le bas).
  https://www.stashofcode.fr/presence-dun-point-dans-un-secteur-angulaire/#more-328
  M    Point dont la position est à tester (point resultant a tester)
  O    Point sommet du secteur (point generateur)
  A    Point de départ du secteur (point adjacent inferieur)
  B    Point d'arrivée du secteur (point adjacent superieur)
  sortie : Booléen indiquant si le point est dans le secteur ou non.
  Exemple :
  B=point(2.0,1.0)
  O=point(2.5,2.5)
  A=point(5.0,5.0)
  M=point(5.0,4.0)
  inSector(M, O, A, B)
=#

function inSector(M, O, A, B)

    cpAB = (A.y - O.y) * (B.x - O.x) - (A.x - O.x) * (B.y - O.y)
    cpAM = (A.y - O.y) * (M.x - O.x) - (A.x - O.x) * (M.y - O.y)
    cpBM = (B.y - O.y) * (M.x - O.x) - (B.x - O.x) * (M.y - O.y)

    if (cpAB > 0)
        if ((cpAM > 0) && (cpBM < 0))
            return true
        else
            return false
        end
    else
        if (!((cpAM < 0) && (cpBM > 0)))
            return true
        else
            return false
        end
    end
end

function inCone(pOrg, pDeb, pFin, pCur)
    # pOrg : point origine du cone (la ou il est pointe)
    # pDeb : point depart du cone (point du rayon [pOrg,pDeb])
    # pFin : point final du cone (point du rayon [pOrg,pFin])
    # pCur : point courant a tester
    # retourne VRAI si pCur est dans le cone pDeb-pFin-pOrg, FAUX sinon

    cp_pDeb_pFin = (pDeb.x - pOrg.x) * (pFin.y - pOrg.y) - (pDeb.y - pOrg.y) * (pFin.x - pOrg.x)
    cp_pDeb_pCur = (pDeb.x - pOrg.x) * (pCur.y - pOrg.y) - (pDeb.y - pOrg.y) * (pCur.x - pOrg.x)
    cp_pFin_pCur = (pFin.x - pOrg.x) * (pCur.y - pOrg.y) - (pFin.y - pOrg.y) * (pCur.x - pOrg.x)

    if (cp_pDeb_pFin > 0)
        if ((cp_pDeb_pCur >= 0) && (cp_pFin_pCur <= 0))
            return true
        else
            return false
        end
    else
        if (!((cp_pDeb_pCur < 0) && (cp_pFin_pCur > 0)))
            return true
        else
            return false
        end
    end
end

function inCone1VersZ(pOrg, pDeb, pFin, pCur)
    return inCone(pOrg, pDeb, pFin, pCur)
end

function inCone2Vers0(pOrg, pDeb, pFin, pCur)
    return !inCone(pOrg, pDeb, pFin, pCur)
end

# ==============================================================================
# Selectionne les points pour le cone pointe sur le generateur k (pCour) et ouvert vers Y
function selectionPoints(vg, k)
    nbgen = size(vg,1)
    if k==1
        # premier generateur (point predecesseur fictif)
        pPrec = tPoint(vg[k].sRel.y[1], vg[k].sRel.y[2]+1.0)
        pCour = tPoint(vg[k].sRel.y[1], vg[k].sRel.y[2])
        pSuiv = tPoint(vg[k+1].sRel.y[1], vg[k+1].sRel.y[2])
    elseif k==nbgen
        # dernier generateur (point suivant fictif)
        pPrec = tPoint(vg[k-1].sRel.y[1], vg[k-1].sRel.y[2])
        pCour = tPoint(vg[k].sRel.y[1], vg[k].sRel.y[2])
        pSuiv = tPoint(vg[k].sRel.y[1]+1.0, vg[k].sRel.y[2])
    else
        # generateur non extreme
        pPrec = tPoint(vg[k-1].sRel.y[1], vg[k-1].sRel.y[2])
        pCour = tPoint(vg[k].sRel.y[1], vg[k].sRel.y[2])
        pSuiv = tPoint(vg[k+1].sRel.y[1], vg[k+1].sRel.y[2])
    end
#    print("Coordonnees du cone 1 : ")
#    @show pPrec, pCour, pSuiv
    return pPrec, pCour, pSuiv
end

# ==============================================================================
# arrondi la solution correspondant au generateur (pas d'historique donc)
# version avec cone inferieur seulement
function roundingSolution!(vg,k,c1,c2)

    nbvar = length(vg[k].sInt.x)
    nbgen = size(vg,1)

    # --------------------------------------------------------------------------
    # identifie les variables fractionnaires et marquage par valeur sentinelle -1
    for i in 1:nbvar
        if  isapprox(vg[k].sPrj.x[i] , 0.0, atol=1e-3)
            vg[k].sInt.x[i] = 0
        elseif isapprox(vg[k].sPrj.x[i] , 1.0, atol=1e-3)
            vg[k].sInt.x[i] = 1
        else
            vg[k].sInt.x[i] = -1
        end
    end

    # --------------------------------------------------------------------------
    # Selectionne les points pour le cone pointe sur le generateur k (pCour) et ouvert vers Y
    pPrec, pCour, pSuiv = selectionPoints(vg, k)

    # --------------------------------------------------------------------------
    # elabore le point pointeur du cone ouvert vers L
    pC = elaborePointConeOuvertversL(vg, k, pPrec, pSuiv)

    # --------------------------------------------------------------------------
    # Arrondi les valeurs non-entieres d'une solution fractionnaire
    # 1) applique un arrondi qui fixe 1 variable (la premiere venue) par iteration

    vg[k].sInt.y[1] = 0
    vg[k].sInt.y[2] = 0

    # calcule en differentiel par rapport a l'image de la solution fractionnaire
    z1 = vg[k].sPrj.y[1]
    z2 = vg[k].sPrj.y[2]
    pM = tPoint( z1 , z2 )

    print("Depart a arrondir : ")
    @show pPrec, pCour, pSuiv, pM
    if length(vg[k].sPrj.x) ≤ 20 @show vg[k].sPrj.x end
    if length(vg[k].sInt.x) ≤ 20 @show vg[k].sInt.x end

    nbVarNonEntiere = 0
    for i in 1:nbvar
        if vg[k].sInt.x[i] == -1

            # la variable i est non entiere
            nbVarNonEntiere += 1

            # defalque la contribution fractionnaire de la variable aux objectifs => pose x[i]=0
            pM = tPoint( z1 - vg[k].sPrj.x[i] * c1[i] , z2 - vg[k].sPrj.x[i] * c2[i])

            if inSector(pM, pCour, pSuiv, pPrec)
                # le point pM obtenu est hors du cone => valide x[i]=1
                vg[k].sInt.x[i] = 1
                z1 = z1 - (vg[k].sPrj.x[i] * c1[i]) + c1[i]
                z2 = z2 - (vg[k].sPrj.x[i] * c2[i]) + c2[i]
                @printf("  x[%2d]=1 |  : [ %12.5f , %12.5f ] \n",i, z1, z2)
            else
                # le point pM obtenu est dans le cone => valide x[i]=0
                vg[k].sInt.x[i] = 0
                z1 = z1 - (vg[k].sPrj.x[i] * c1[i])
                z2 = z2 - (vg[k].sPrj.x[i] * c2[i])
                @printf("  x[%2d]=0 |  : [ %12.5f , %12.5f ] \n",i, z1, z2)
            end
        end

        # calcule la performance de la solution entiere
        vg[k].sInt.y[1] += vg[k].sInt.x[i] * c1[i]
        vg[k].sInt.y[2] += vg[k].sInt.x[i] * c2[i]
    end

    if length(vg[k].sInt.x) ≤ 20 @show vg[k].sInt.x end
    @printf("→ #round : %4d → [ %5d , %5d ] ", nbVarNonEntiere, vg[k].sInt.y[1], vg[k].sInt.y[2])
    push!(XInt,vg[k].sInt.y[1])
    push!(YInt,vg[k].sInt.y[2])

end

# ==============================================================================
# arrondi la solution correspondant au generateur (pas d'historique donc)
# version avec cone inferieur et superieur
function roundingSolutionnew24!(vg,k,c1,c2)

    nbvar = length(vg[k].sInt.x)
    nbgen = size(vg,1)

    # --------------------------------------------------------------------------
    # identifie les variables fractionnaires et marquage par valeur sentinelle -1
    for i in 1:nbvar
        if  isapprox(vg[k].sPrj.x[i] , 0.0, atol=1e-3)
            vg[k].sInt.x[i] = 0
        elseif isapprox(vg[k].sPrj.x[i] , 1.0, atol=1e-3)
            vg[k].sInt.x[i] = 1
        else
            vg[k].sInt.x[i] = -1
        end
    end

    # --------------------------------------------------------------------------
    # Selectionne les points pour le cone pointe sur le generateur k (pCour) et ouvert vers Y
    pPrec, pCour, pSuiv = selectionPoints(vg, k)

    # --------------------------------------------------------------------------
    # elabore le point pointeur du cone ouvert vers L
    pC = elaborePointConeOuvertversL(vg, k, pPrec, pSuiv)

    # --------------------------------------------------------------------------
    # Arrondi les valeurs non-entieres d'une solution fractionnaire
    # 1) applique un arrondi qui fixe 1 variable (la premiere venue) par iteration

    vg[k].sInt.y[1] = 0
    vg[k].sInt.y[2] = 0

    # calcule en differentiel par rapport a l'image de la solution fractionnaire
    z1 = vg[k].sPrj.y[1]
    z2 = vg[k].sPrj.y[2]
    pM = tPoint( z1 , z2 )

#    print("Depart a arrondir : ")
#    @show pPrec, pCour, pSuiv, pM
#    if length(vg[k].sPrj.x) ≤ 20 @show vg[k].sPrj.x end
#    if length(vg[k].sInt.x) ≤ 20 @show vg[k].sInt.x end

    verbose ? @printf("  %2dR : [ %8.2f , %8.2f ] ", k, z1, z2) : nothing

    nbVarNonEntiere = 0
    for i in 1:nbvar
        if vg[k].sInt.x[i] == -1

            # la variable i est non entiere
            nbVarNonEntiere += 1

            # defalque la contribution fractionnaire de la variable aux objectifs => pose x[i]=0
            pM0 = tPoint( z1 - vg[k].sPrj.x[i] * c1[i] , z2 - vg[k].sPrj.x[i] * c2[i])
            pM1 = tPoint( z1 - vg[k].sPrj.x[i] * c1[i]+c1[i] , z2 - vg[k].sPrj.x[i] * c2[i]+c2[i])

            if !inCone1VersZ(pCour, pSuiv, pPrec, pM0)
                # le point pM0 obtenu est hors du cone inferieur => valide x[i]=1
                vg[k].sInt.x[i] = 1
                z1 = z1 - (vg[k].sPrj.x[i] * c1[i]) + c1[i]
                z2 = z2 - (vg[k].sPrj.x[i] * c2[i]) + c2[i]
#                @printf("  x[%2d]=1 |  : [ %12.5f , %12.5f ] \n",i, z1, z2)
            elseif inCone2Vers0(pC, pSuiv, pPrec, pM0)
                # le point pM0 obtenu est dans le cone => valide x[i]=0
                vg[k].sInt.x[i] = 0
                z1 = z1 - (vg[k].sPrj.x[i] * c1[i])
                z2 = z2 - (vg[k].sPrj.x[i] * c2[i])
#                @printf("  x[%2d]=0 |  : [ %12.5f , %12.5f ] \n",i, z1, z2)
            elseif inCone2Vers0(pC, pSuiv, pPrec, pM1)
                # le point pM1 obtenu est dans le cone => valide x[i]=1
                vg[k].sInt.x[i] = 1
                z1 = z1 - (vg[k].sPrj.x[i] * c1[i]) + c1[i]
                z2 = z2 - (vg[k].sPrj.x[i] * c2[i]) + c2[i]
#                @printf("  x[%2d]=0 |  : [ %12.5f , %12.5f ] \n",i, z1, z2)
            else
#                println("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
                # choisir entre pM0 et pM1 selon la distance a la droite [O;OO]
                O = tPoint(vg[k].sRel.y[1], vg[k].sRel.y[2]) # generateur
                OO = tPoint(vg[k].sPrj.y[1], vg[k].sPrj.y[2]) # projection
                # les points pM0 et pM1 obtenus sont hors du cone => choix de 0 ou 1 avec la plus courte distance à O et OO
                distM0 = sqrt((O.x-pM0.x)^2+(O.y-pM0.y)^2)+sqrt((OO.x-pM0.x)^2+(OO.y-pM0.y)^2)
                distM1 = sqrt((O.x-pM1.x)^2+(O.y-pM1.y)^2)+sqrt((OO.x-pM1.x)^2+(OO.y-pM1.y)^2)
                if distM0 < distM1
                    vg[k].sInt.x[i] = 0
                    z1 = z1 - (vg[k].sPrj.x[i] * c1[i])
                    z2 = z2 - (vg[k].sPrj.x[i] * c2[i])
                else
                    vg[k].sInt.x[i] = 1
                    z1 = z1 - (vg[k].sPrj.x[i] * c1[i]) + c1[i]
                    z2 = z2 - (vg[k].sPrj.x[i] * c2[i]) + c2[i]
                end
            end
        end

        # calcule la performance de la solution entiere
        vg[k].sInt.y[1] += vg[k].sInt.x[i] * c1[i]
        vg[k].sInt.y[2] += vg[k].sInt.x[i] * c2[i]
    end

    if length(vg[k].sInt.x) ≤ 20 @show vg[k].sInt.x end
    @printf("→ #round : %4d → [ %5d , %5d ] ", nbVarNonEntiere, vg[k].sInt.y[1], vg[k].sInt.y[2])
    push!(XInt,vg[k].sInt.y[1])
    push!(YInt,vg[k].sInt.y[2])

end

# ==============================================================================
# arrondi la solution correspondant au generateur (pas d'historique donc)
# version avec voisinage et selection d'un voisin selon distance L1 avec generateur
function roundingSolutionNew23!(vg,k,c1,c2)

    nbvar = length(vg[k].sInt.x)
    nbgen = size(vg,1)
    lstIdxFrac =(Int64)[]

    vg[k].sInt.y[1] = 0
    vg[k].sInt.y[2] = 0

    # --------------------------------------------------------------------------
    # identifie les variables fractionnaires et marquage par valeur sentinelle -1
    for i in 1:nbvar
        if  isapprox(vg[k].sPrj.x[i] , 0.0, atol=1e-3)
            vg[k].sInt.x[i] = 0
        elseif isapprox(vg[k].sPrj.x[i] , 1.0, atol=1e-3)
            vg[k].sInt.x[i] = 1
            vg[k].sInt.y[1] += c1[i] # Comptabilise la contribution des var a 1
            vg[k].sInt.y[2] += c2[i]
        else
            vg[k].sInt.x[i] = -1
            push!(lstIdxFrac,i)
        end
    end

    # --------------------------------------------------------------------------
    # Selectionne les points pour le cone pointe sur le generateur k (pCour) et ouvert vers Y
    pPrec, pCour, pSuiv = selectionPoints(vg, k)

    # --------------------------------------------------------------------------
    # elabore le point pointeur du cone ouvert vers L
    pC = elaborePointConeOuvertversL(vg, k, pPrec, pSuiv)

    # --------------------------------------------------------------------------
    # Arrondi les valeurs non-entieres d'une solution fractionnaire
    # 1) applique un arrondi qui fixe 1 variable (la premiere venue) par iteration

    # calcule en differentiel par rapport a l'image de la solution fractionnaire
    z1 = vg[k].sPrj.y[1]
    z2 = vg[k].sPrj.y[2]
    pM = tPoint( z1 , z2 )

#    if length(vg[k].sPrj.x) ≤ 20 @show vg[k].sPrj.x end
#    if length(vg[k].sInt.x) ≤ 20 @show vg[k].sInt.x end

    nbVarNonEntiere = length(lstIdxFrac)
    # traite toutes les variables non entieres

#=
    @printf("\n\n")
    @show c1
    @show c2
    @printf("VARIABLES FIXEES : \n")
    @show vg[k].sPrj.x
    @show vg[k].sInt.x
    @show vg[k].sInt.y
    @show pM
    @show lstIdxFrac
=#
    verbose ? @printf("  %2dR : [ %8.2f , %8.2f ] ", k, z1, z2) : nothing

    # traitement de toutes les variables fractionnaires une a une
    for i= 1:nbVarNonEntiere
#        @printf("\n")

        distMin = Inf; ind = -1; val = -1
        # voisinage a 1 variable sur toutes les variables
        for iFrac in [i for i=1:nbvar if vg[k].sInt.x[i] == -1]

            # Evalue x[iFrac]=0
            # defalque la contribution fractionnaire de la variable aux objectifs => pose x[iFrac]=0
            pM = tPoint( z1 - c1[iFrac] * vg[k].sPrj.x[iFrac] , z2 - c2[iFrac] * vg[k].sPrj.x[iFrac])
            if inCone1VersZ(pCour, pSuiv, pPrec, pM)
                distL1 = abs(vg[k].sRel.y[1]-pM.x) + abs(vg[k].sRel.y[2]-pM.y)
                if distL1 < distMin
                    distMin = distL1; ind = iFrac; val = 0
                end
            end

            # Evalue x[iFrac]=1
            # ajoute la contribution entiere de la variable aux objectifs => pose x[iFrac]=1
            pM = tPoint( z1 - c1[iFrac] * vg[k].sPrj.x[iFrac] + c1[iFrac], z2 - c2[iFrac] * vg[k].sPrj.x[iFrac] + c2[iFrac])
            distL1 = abs(vg[k].sRel.y[1]-pM.x) + abs(vg[k].sRel.y[2]-pM.y)
            if distL1 < distMin
                distMin = distL1; ind = iFrac; val = 1
            end
        end

        # applique le voisin le plus interessant trouve
        vg[k].sInt.x[ind] = val
        vg[k].sInt.y[1] += vg[k].sInt.x[ind] * c1[ind] * val
        vg[k].sInt.y[2] += vg[k].sInt.x[ind] * c2[ind] * val
        z1 = z1 - (vg[k].sPrj.x[ind] * c1[ind]) + c1[ind] * val
        z2 = z2 - (vg[k].sPrj.x[ind] * c2[ind]) + c2[ind] * val
    #    @printf("  x[%2d]=%2d |  : [ %12.5f , %12.5f ] \n",ind, val, z1, z2)

    end

#    if length(vg[k].sInt.x) ≤ 20 @show vg[k].sInt.x end
    @printf("→ #round : %4d → [ %5d , %5d ] ", nbVarNonEntiere, vg[k].sInt.y[1], vg[k].sInt.y[2])
    push!(XInt,vg[k].sInt.y[1])
    push!(YInt,vg[k].sInt.y[2])

end

# ==============================================================================
# projecte la solution entiere correspondant au generateur k et test d'admissibilite
function projectingSolution!(vg,k,A,c1,c2,λ1,λ2)

    # --------------------------------------------------------------------------
    # Projete la solution entiere sur le polytope X avec norme-L1

#    fPrj, vg[k].sPrj.x = Δ2SPA(A,vg[k].sInt.x)
    fPrj, vg[k].sPrj.x = Δ2SPAbis(A,vg[k].sInt.x, c1, c2,k,λ1,λ2)

    # Nettoyage de la valeur de vg[k].sPrj.x et calcul du point bi-objectif
    # reconditionne les valeurs 0 et 1 et arrondi les autres valeurs
    nettoyageSolution!(vg[k].sPrj.x)
#    verbose ? @printf("  %2dP : fPrj = %8.2f  ",k, round(fPrj, digits=2)) : nothing

    # recalcule la solution au regard des 2 objectifs
    vg[k].sPrj.y[1], vg[k].sPrj.y[2] = evaluerSolution(vg[k].sPrj.x, c1, c2)
    verbose ? @printf("  %2dP : [ %8.2f , %8.2f ] ",k, vg[k].sPrj.y[1], vg[k].sPrj.y[2]) : nothing

    push!(XProj, vg[k].sPrj.y[1])
    push!(YProj, vg[k].sPrj.y[2])

    # ----------------------------------------------------------------
    # Teste si la projection est admissible

    if estAdmissible(vg[k].sPrj.x)
        vg[k].sInt.x = deepcopy(vg[k].sPrj.x)
        vg[k].sInt.y[1] = vg[k].sPrj.y[1]
        vg[k].sInt.y[2] = vg[k].sPrj.y[2]
        vg[k].sFea = true
        push!(XFeas, vg[k].sPrj.y[1])
        push!(YFeas, vg[k].sPrj.y[2])
        @printf("→ Admissible "); print("                       ")
    else
        vg[k].sFea = false
        @printf("→ x          "); print("                       ")
        # prepare pour l'iteration suivante
#        vg[k].xRlx = deepcopy(vg[k].sPrj.x) !!!!!!!!!!!!!
    end

end

# ==============================================================================
# Calcule la direction d'interet du nadir vers le milieu de segment reliant deux points generateurs
function calculerDirections(L,vg)

    nbgen = size(vg,1)
    for k in 2:nbgen

        n1 = L[end].y[1]
        n2 = L[1].y[2]

        x1,y1 = vg[k-1].sRel.y[1], vg[k-1].sRel.y[2]
        x2,y2 = vg[k].sRel.y[1], vg[k].sRel.y[2]
        xm=(x1+x2)/2.0
        ym=(y1+y2)/2.0
        Δx = abs(n1-xm)
        Δy = abs(n2-ym)
        λ1 =  1 - Δx / (Δx+Δy)
        λ2 =  1 - Δy / (Δx+Δy)
        @printf("  x1= %7.2f   y1= %7.2f \n",x1,y1)
        @printf("  x2= %7.2f   y2= %7.2f \n",x2,y2)
        @printf("  Δx= %7.2f    Δy= %7.2f \n",Δx,Δy)
        @printf("  λ1= %6.5f    λ2= %6.5f \n",λ1,λ2)
        plot(n1, n2, xm, ym, linestyle="-", color="blue", marker="+")
        annotate("",
                 xy=[xm;ym],# Arrow tip
                 xytext=[n1;n2], # Text offset from tip
                 arrowprops=Dict("arrowstyle"=>"->"))
        println("")
    end

end

# ==============================================================================
# Calcule la direction d'interet du nadir vers un point generateur
function calculerDirections2(L,vg)

    nbgen = size(vg,1)
    λ1=Vector{Float64}(undef, nbgen)
    λ2=Vector{Float64}(undef, nbgen)
    for k in 1:nbgen

        n1 = L[end].y[1]
        n2 = L[1].y[2]

        xm=vg[k].sRel.y[1]
        ym=vg[k].sRel.y[2]
        Δx = abs(n1-xm)
        Δy = abs(n2-ym)
        λ1[k] =  1 - Δx / (Δx+Δy)
        λ2[k] =  1 - Δy / (Δx+Δy)
        @printf("  k= %3d   ",k)
        @printf("  xm= %7.2f   ym= %7.2f ",xm,ym)
        @printf("  Δx= %8.2f    Δy= %8.2f ",Δx,Δy)
        @printf("  λ1= %6.5f    λ2= %6.5f \n",λ1[k],λ2[k])
        plot(n1, n2, xm, ym, linestyle="-", color="blue", marker="+")
        annotate("",
                 xy=[xm;ym],# Arrow tip
                 xytext=[n1;n2], # Text offset from tip
                 arrowprops=Dict("arrowstyle"=>"->"))
        #println("")
    end
    return λ1, λ2
end

# ==============================================================================
# applique une perturbation sur la solution entiere faisant l'objet d'un cycle
function perturbSolution!(vg,k,c1,c2)
    # nombre de variables maximum a considerer
    T = 20
    # nombre effectif de variables a flipper
    TT = rand(T/2:3*T/2)

    # liste des candidats (valeur, indice) et tri decroissant
    nbvar = length(vg[k].sInt.x)
    candidats=[( abs( vg[k].sPrj.x[i] - vg[k].sInt.x[i] ) , i ) for i=1:nbvar]
    sort!(candidats, rev=true, by = x -> x[1])
#    sort!(candidats,  by = x -> x[1])

#@show vg[k].sPrj.x
#@show vg[k].sInt.x
#@show candidats
#@show TT
#@show nbvar

    i = 1
    while (i<= nbvar) && (i<=TT)
        j=candidats[i][2]
#        @show candidats[i][2]
        if vg[k].sInt.x[j] == 0
            vg[k].sInt.x[j] = 1
            vg[k].sInt.y[1] = vg[k].sInt.y[1] + c1[j]
            vg[k].sInt.y[2] = vg[k].sInt.y[2] + c2[j]
        else
            vg[k].sInt.x[j] = 0
            vg[k].sInt.y[1] = vg[k].sInt.y[1] - c1[j]
            vg[k].sInt.y[2] = vg[k].sInt.y[2] - c2[j]
        end
        i+=1
    end
    @printf("  %2dC : [ %8.2f , %8.2f ] \n", k, vg[k].sInt.y[1], vg[k].sInt.y[2])
    push!(XPert,vg[k].sInt.y[1]); push!(YPert,vg[k].sInt.y[2])
#    @show vg[k].sInt.x
    return nothing
end

# ==============================================================================
# applique une perturbation sur la solution entiere faisant l'objet d'un cycle
function perturbSolution30!(vg,k,c1,c2)

    # liste des candidats (valeur, indice) et tri decroissant
    nbvar = length(vg[k].sInt.x)
    idxTilde0, idxTilde1 = split01(vg[k].sInt.x)

#    candidats=[( abs( vg[k].sPrj.x[i] - vg[k].sInt.x[i] ) , i ) for i=1:nbvar if vg[k].sPrj.x[i]>0 && vg[k].sPrj.x[i]<1]
    candidats=[( vg[k].sPrj.x[i] , i ) for i=1:nbvar if vg[k].sPrj.x[i]>0 && vg[k].sPrj.x[i]<1]
#    sort!(candidats, rev=true, by = x -> x[1])

#@show vg[k].sPrj.x
#@show vg[k].sInt.x
#@show candidats
#@show nbvar

seq = randperm(length(candidats)) # melange les candidats afin d'avoir une composante variee
etat = vg[k].sInt.x[ candidats[seq[1]][2] ] # etat 0 ou 1 de la premiere variable candidate
for i = 1:length(candidats)
    j=candidats[seq[i]][2]
    #@show candidats[seq[i]][2]
    if etat == 0
        if vg[k].sInt.x[j] == 0
            vg[k].sInt.x[j] = 1
            vg[k].sInt.y[1] = vg[k].sInt.y[1] + c1[j]
            vg[k].sInt.y[2] = vg[k].sInt.y[2] + c2[j]
        end
    else
        if vg[k].sInt.x[j] == 1
            vg[k].sInt.x[j] = 0
            vg[k].sInt.y[1] = vg[k].sInt.y[1] - c1[j]
            vg[k].sInt.y[2] = vg[k].sInt.y[2] - c2[j]
        end
    end
    etat=(etat+1)%2
end

    @printf("  %2dC : [ %8.2f , %8.2f ] \n", k, vg[k].sInt.y[1], vg[k].sInt.y[2])
    push!(XPert,vg[k].sInt.y[1]); push!(YPert,vg[k].sInt.y[2])
#    @show vg[k].sInt.x

    return nothing
end

# ==============================================================================
# point d'entree principal

function GM( fname::String,
             tailleSampling::Int64,
             maxTrial::Int64,
             maxTime::Int64
           )

    @assert tailleSampling>=3 "Erreur : Au moins 3 sont requis"

    @printf("0) instance et parametres \n\n")
    verbose ? println("  instance = $fname | tailleSampling = $tailleSampling | maxTrial = $maxTrial | maxTime = $maxTime\n\n") : nothing

    # chargement de l'instance numerique ---------------------------------------
    c1, c2, A = loadInstance2SPA(fname) # instance numerique de SPA
    nbctr = size(A,1)
    nbvar = size(A,2)
    nbobj = 2

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    @printf("1) calcule les etendues de valeurs sur les 2 objectifs\n\n")

    # calcule la valeur optimale relachee de f1 seule et le point (z1,z2) correspondant
    f1RL, xf1RL = computeLinearRelax2SPA(nbvar, nbctr, A, c1, c2, typemax(Int), 1) # opt fct 1
    minf1RL, maxf2RL = evaluerSolution(xf1RL, c1, c2)

    # calcule la valeur optimale relachee de f2 seule et le point (z1,z2) correspondant
    f2RL, xf2RL = computeLinearRelax2SPA(nbvar, nbctr, A, c1, c2, typemax(Int), 2) # opt fct 2
    maxf1RL, minf2RL = evaluerSolution(xf2RL, c1, c2)

    verbose ? @printf("  f1_min=%8.2f ↔ f1_max=%8.2f (Δ=%.2f) \n",minf1RL, maxf1RL, maxf1RL-minf1RL) : nothing
    verbose ? @printf("  f2_min=%8.2f ↔ f2_max=%8.2f (Δ=%.2f) \n\n",minf2RL, maxf2RL, maxf2RL-minf2RL) : nothing

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    @printf("2) calcule les generateurs par e-contrainte alternant minimiser z1 et z2\n\n")

    nbgen, L = calculGenerateurs(A, c1, c2, tailleSampling, minf1RL, maxf2RL, maxf1RL, minf2RL)

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    # allocation de memoire pour la structure de donnees -----------------------

    vg = allocateDatastructure(nbgen, nbvar, nbobj)

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    @printf("3) place L dans structure et verifie l'admissibilite de chaque generateur\n\n")

    for k=1:nbgen

        verbose ? @printf("  %2d  : [ %8.2f , %8.2f ] ", k, L[k].y[1], L[k].y[2]) : nothing

        # copie de l'ensemble bornant inferieur dans la stru de donnees iterative ---
        ajouterX0!(vg, k, L[k])

        # test d'admissibilite et marquage de la solution le cas echeant -------
        if estAdmissible(vg[k].sRel.x)
            ajouterXtilde!(vg, k, convert.(Int, vg[k].sRel.x), convert.(Int, L[k].y))
            vg[k].sFea   = true
            push!(XFeas,vg[k].sInt.y[1])
            push!(YFeas,vg[k].sInt.y[2])
            verbose ? @printf("→ Admissible \n") : nothing
        else
            vg[k].sFea   = false
            verbose ? @printf("→ x          \n") : nothing
        end

    end
    verbose ? println("") : nothing

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    # Sortie graphique

    figure("Gravity Machine",figsize=(6.5,5))
    #xlim(25000,45000)
    #ylim(20000,40000)
    xlabel(L"z^1(x)")
    ylabel(L"z^2(x)")
    PyPlot.title("Cone | 1 rounding | 2-$fname")

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    # calcule les directions (λ1,λ2) pour chaque generateur a utiliser lors des projections
    λ1,λ2 = calculerDirections2(L,vg)

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

    @printf("4) terraformation generateur par generateur \n\n")

    for k in [i for i in 1:nbgen if !isFeasible(vg,i)]
        temps = time()
        trial = 0
        H =(Vector{Int64})[]

#perturbSolution30!(vg,k,c1,c2)

        # rounding solution : met a jour sInt dans vg --------------------------
        #roundingSolution!(vg,k,c1,c2)  # un cone
        #roundingSolutionnew24!(vg,k,c1,c2) # deux cones
        roundingSolutionNew23!(vg,k,c1,c2) # un cone et LS sur generateur

        push!(H,[vg[k].sInt.y[1],vg[k].sInt.y[2]])
        println("   t=",trial,"  |  Tps=", round(time()- temps, digits=4))

        while !(t1=isFeasible(vg,k)) && !(t2=isFinished(trial, maxTrial)) && !(t3=isTimeout(temps, maxTime))

            trial+=1

            # projecting solution : met a jour sPrj, sInt, sFea dans vg --------
            projectingSolution!(vg,k,A,c1,c2,λ1,λ2)
            println("   t=",trial,"  |  Tps=", round(time()- temps, digits=4))

            if !isFeasible(vg,k)

                # rounding solution : met a jour sInt dans vg --------------------------
                #roundingSolution!(vg,k,c1,c2)
                #roundingSolutionnew24!(vg,k,c1,c2)
                roundingSolutionNew23!(vg,k,c1,c2)
                println("   t=",trial,"  |  Tps=", round(time()- temps, digits=4))

                # test detection cycle sur solutions entieres ------------------
                cycle = [vg[k].sInt.y[1],vg[k].sInt.y[2]] in H
                if (cycle == true)
                    println("CYCLE!!!!!!!!!!!!!!!")
                    # perturb solution
                    perturbSolution30!(vg,k,c1,c2)
                end
                push!(H,[vg[k].sInt.y[1],vg[k].sInt.y[2]])

            end
        end
        if t1
            println("   feasible \n")
        elseif t2
            println("   maxTrial \n")
        elseif t3
            println("   maxTime \n")
        end

    end

    println("");

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

    @printf("5) Extraction des resultats\n\n")

    for k=1:nbgen
        verbose ? @printf("  %2d  : [ %8.2f , %8.2f ] ", k, vg[k].sInt.y[1],vg[k].sInt.y[2]) : nothing
        # test d'admissibilite et marquage de la solution le cas echeant -------
        if vg[k].sFea
            verbose ? @printf("→ Admissible \n") : nothing
        else
            verbose ? @printf("→ x          \n") : nothing
        end
    end

    # allocation de memoire pour les ensembles bornants ------------------------
    U = Vector{tSolution{Int64}}(undef,nbgen)
    for j = 1:nbgen
        U[j] = tSolution{Int64}(zeros(Int64,nbvar),zeros(Int64,nbobj))
    end

    # ==========================================================================
    @printf("6) Edition des resultats \n\n")

#    figure("Gravity Machine",figsize=(6.5,5))
    #xlim(25000,45000)
    #ylim(20000,40000)
#    xlabel(L"z^1(x)")
#    ylabel(L"z^2(x)")
    # Donne les points relaches initiaux ---------------------------------------
#    scatter(xLf1,yLf1,color="blue", marker="x")
#    scatter(xLf2,yLf2,color="red", marker="+")
    graphic ? scatter(xL,yL,color="blue", marker="x", label = L"y \in L") : nothing

    # Donne les points entiers -------------------------------------------------
    graphic ? scatter(XInt,YInt,color="orange", marker="s", label = L"y"*" rounded") : nothing
#    @show XInt
#    @show YInt

    # Donne les points apres projection Δ(x,x̃) ---------------------------------
    graphic ? scatter(XProj,YProj, color="red", marker="x", label = L"y"*" projected") : nothing
#    @show XProj
#    @show YProj

    # Donne les points admissibles ---------------------------------------------
    graphic ? scatter(XFeas,YFeas, color="green", marker="o", label = L"y \in F") : nothing
#    @show XFeas
#    @show YFeas

    # Donne l'ensemble bornant primal obtenu + la frontiere correspondante -----
    X_EBP_frontiere, Y_EBP_frontiere, X_EBP, Y_EBP = ExtractEBP(XFeas, YFeas)
    plot(X_EBP_frontiere, Y_EBP_frontiere, color="green", markersize=3.0, marker="x")
    scatter(X_EBP, Y_EBP, color="green", s = 150, alpha = 0.3, label = L"y \in U")
    @show X_EBP
    @show Y_EBP

    # Donne les points qui ont fait l'objet d'une perturbation -----------------
     scatter(XPert,YPert, color="magenta", marker="s", label ="pertub")

    # Donne les points non-domines exacts de cette instance --------------------
     XN,YN = loadNDPoints2SPA(fname)
     plot(XN, YN, color="black", linewidth=0.75, marker="+", markersize=1.0, linestyle=":", label = L"y \in Y_N")
     scatter(XN, YN, color="black", marker="+")
#    @show X_Y_N
#    @show Y_Y_N

    # Affiche le cadre avec les legendes des differents traces -----------------
    legend(bbox_to_anchor=[1,1], loc=0, borderaxespad=0, fontsize = "x-small")
    #PyPlot.title("Cone | 1 rounding | 2-$fname")

end

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

#testidee()
#@time GM("sppaa02.txt", 6, 20, 20)
#@time GM("sppnw03.txt", 6, 20, 20)
#@time GM("sppnw10.txt", 6, 20, 20)
@time GM("didactic5.txt", 5, 5, 10)
#@time GM("sppnw29.txt", 6, 30, 20)
nothing

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine