text:
# # the ideas are based on https://github.com/tkf/Bifurcations.jl/blob/master/src/diffeq.jl
# # I would just grab DiffEqDiffTools.uJacobianWrapper
# # https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/1c86896dfcc0c383ef501014c56a3f8dfc8a63bb/src/nlsolve/newton.jl
# using DiffEqBase: AbstractODEProblem
# using DiffEqOperators
# using Setfield
# using OrdinaryDiffEq
#
# const DEP{iip} = AbstractODEProblem{uType, tType, iip} where {uType, tType}
#
# struct DiffEqWrapper{P, L}
# 	de_prob::P
# 	param_axis::L
# end
#
# mutable struct JacWrap{Tj, Tp, Tx}
# 	J::Tj	# jacobian
# 	p::Tp	# parameter
# 	#TODO: remove this
# 	x::Tx	# vecteur where the jacobian is evaluated
# end
#
# # callable function to compute the Jacobian J(x) = dF(x)
# function (Jac::JacWrap)(x)
# 	if Jac.J isa DiffEqBase.AbstractDiffEqLinearOperator
# 		Jac.x .= x
# 		update_coefficients(Jac.J, x, Jac.p, zero(eltype(x)))
# 		# this is a Hack because update_coefficients does not do it
# 		Jac.J.u .= x
# 		norm(Jac.J.u - x) ≈ 0 ? printstyled(color=:green,"--> Jac Updated\n") : printstyled(color=:red,"--> Jac Not updated!!, $(norm(Jac.J.u - x))\n")
# 	else
# 		@error "Not a Matrix Free Operator!"
# 	end
# 	# needs to return the jacobian so we can call linsolve on it
# 	Jac.J
# end
#
# # callable function to compute the Jacobian J(x, p)
# function (Jac::JacWrap{Tj, Tp, Tx})(x, p::Tp) where {Tj, Tp, Tx}
# 	Jac.p = p
# 	Jac(x)
# end
#
# ####################################################################################################
# # we define a linearsolver to use our interface based on the one in DiffEqBase in linear_nonlinear.jl
#
# struct DiffEqWrapLS <: AbstractLinearSolver
# 	linsolve
# end
#
# function (l::DiffEqWrapLS)(J, rhs)
# 	out = similar(rhs)
#
# 	# tol is required for DefaultLinSolve() to work
# 	l.linsolve(out, J, rhs; tol = 1e-4, verbose = true)
#
# 	return out, true, 1
# end
# ####################################################################################################
# # Nonlinear solvers
#
# function ContinuationProblem(
# 		deprob::DEP{iip}, param_axis::Lens;
# 		kwargs0...) where iip
# 	de_prob = deepcopy(deprob)
# 	x0 = copy(de_prob.u0)
# 	p0 = de_prob.p
# 	@assert !(typeof(x0) <: Number) "We need array like structures for the state space, consider using [u0]"
# 	pbwrap = DiffEqWrapper(de_prob, param_axis)
#
# 	# we extract the vector field
# 	if de_prob isa DEP{false}
# 		fnewton = x -> de_prob.f(x, p0, de_prob.tspan[1])
# 		f = (x, p) -> de_prob.f(x, set(p0, param_axis, p), de_prob.tspan[1])
# 	else
# 		fnewton = x -> (out = similar(x);de_prob.f(out, x, p0, de_prob.tspan[1]);out)
# 		f =  (x, p) -> (out = similar(x);de_prob.f(out, x, set(p0, param_axis, p), de_prob.tspan[1]);out)
# 	end
#
# 	# we extract the jacobian (operator)
# 	if de_prob.f.jac === nothing
# 		@warn "No jacobian is provided by ODEProblem"
# 		dfnewton = nothing
# 		df = nothing
# 	else
# 		@warn "jacobian is provided by ODEProblem"
# 		J, W = OrdinaryDiffEq.build_J_W(ImplicitEuler(), x0, x0, p0, 1., 0., de_prob.f, eltype(x0), Val{}(iip))
# 		@show typeof(J), typeof(W)
# 		dfnewton = J
# 		df = nothing
# 	end
# 	return (fnewton = fnewton, dfnewton = dfnewton, f = f, df = df, x0 = x0)
# end
#
#
# function newton(de_prob::DEP, options:: NewtonPar{T}; normN = norm, linsolver = DiffEqBase.DefaultLinSolve()) where {T}
# 	@unpack fnewton, dfnewton, x0 = ContinuationProblem(de_prob, @lens _.p)
# 	if dfnewton == nothing
# 		@warn "newton without Jacobian"
# 		return newton(fnewton, x0, options)
# 	else
# 		@warn "newton with Jacobian"
# 		options = @set options.linsolver = DiffEqWrapLS(linsolver)
# 		J = JacWrap(dfnewton, de_prob.p, copy(x0))
# 		return newton(fnewton, J, x0, options)
# 	end
# end
#
# function continuation(de_prob::DEP, param_axis::Lens, contParams::ContinuationPar{T, S, E};kwargs...) where {T, S, E}
# 	@unpack f, df, x0 = ContinuationProblem(de_prob, param_axis)
# 	p = Setfield.get(de_prob.p, param_axis)
# 	continuation(f, x0, p, contParams; kwargs...)
# end


text:
function selected(tree_view::GtkTreeView, list::GtkTreeStore)
    selmodel = Gtk.GAccessor.selection(tree_view)
    if hasselection(selmodel)
        iter = selected(selmodel)
        return list[iter]
    end
    return nothing
end

#select the first entry that is equal to v
function select_value(tree_view::GtkTreeView,list::GtkTreeStore,v)
    selmodel = Gtk.GAccessor.selection(tree_view)
    for i = 1:length(list)
        if list[i] == v
            partialsort!(selmodel, Gtk.iter_from_index(list, i))
            return
        end
    end
end

expand_root(tree_view::GtkTreeView) = Gtk.expand_to_path(tree_view, Gtk.treepath("0"))

gdk_keyval_name(val) = unsafe_string(
    ccall((:gdk_keyval_name,libgtk),Ptr{UInt8},(Cuint,),val),
true)

ext: jl


text:
# This file defines some struct types used in the FrameConnection package.

using DataFrames

"""
    SMD(connectID::Vector{Int}, x::Vector{Float32}, y::Vector{Float32},
        x_se::Vector{Float32}, y_se::Vector{Float32}, framenum::Vector{Int},
        nframes::Int, datasetnum::Vector{Int}, ndatasets::Int)

Single Molecule Data structure with information about localizations.
"""
mutable struct SMD
    connectID::Vector{Int}
    x::Vector{Float32}
    y::Vector{Float32}
    x_se::Vector{Float32}
    y_se::Vector{Float32}
    framenum::Vector{Int}
    nframes::Int
    datasetnum::Vector{Int}
    ndatasets::Int
    SMD() = new()
end
function SMD(nlocalizations::Int)
    smd = SMD()
    smd.x = Vector{Float32}(undef, nlocalizations)
    smd.y = Vector{Float32}(undef, nlocalizations)
    smd.x_se = Vector{Float32}(undef, nlocalizations)
    smd.y_se = Vector{Float32}(undef, nlocalizations)
    smd.framenum = Vector{Int}(undef, nlocalizations)
    smd.datasetnum = Vector{Int}(undef, nlocalizations)

    return smd
end
function SMD(data::DataFrames.DataFrame)
    smd = SMD()
    smd.datasetnum = Int.(data[:, 1])
    smd.framenum = Int.(data[:, 2])
    smd.x = Float32.(data[:, 3])
    smd.y = Float32.(data[:, 4])
    smd.x_se = Float32.(data[:, 5])
    smd.y_se = Float32.(data[:, 6])
    smd.nframes = Int(maximum(smd.framenum))
    smd.ndatasets = Int(length(unique(smd.datasetnum)))

    return smd
end

"""
    ParamStruct(initialdensity::Vector{Float32},
        k_on::Float32, k_off::Float32, k_bleach::Float32, p_miss::Float32,
        nsigmadev::Float32, maxframegap::Int, nnearestclusters::Int)

Structure of parameters needed for frame-connection.

# Fields
-`initialdensity`: Density of emitters at the start of the experiment.
                   (see estimatedensities()) (emitters/pixel^2)
-`nnearestclusters`: Number of nearest preclusters used for local density
                     estimates. (default = 2)(see estimatedensities())
-`k_on`: Rate at which dark emitters convert to the visible state.
         (see estimateparams())(1/frame)
-`k_off`: Rate at which visible emitters are converted to the reversible dark
          state. (see estimateparams())(1/frame)
-`k_bleach`: Rate at which visible emitters are irreversibly photobleached.
             (see estimateparams())(1/frame)
-`p_miss`: Probability of missing a localization of a visible emitter.
-`nsigmadev`: Multiplier of localization errors that defines a pre-clustering
              distance threshold. (default = 5)(see precluster())(pixels)
-`maxframegap`: Maximum frame gap between temporally adjacent localizations in
                a precluster. (default = 5)(see precluster())(frames)
-`nmaxnn`: Maximum number of nearest-neighbors inspected for precluster
           membership.  Ideally, this would be set to inf, but that's not
           feasible for most data. (default = 2)(see precluster())

"""
mutable struct ParamStruct
    initialdensity::Vector{Float32}
    nnearestclusters::Int
    k_on::Float32
    k_off::Float32
    k_bleach::Float32
    p_miss::Float32
    nsigmadev::Float32
    maxframegap::Int
    nmaxnn::Int
end
ParamStruct() = ParamStruct([], 2, 0, 0, 0, 0, 5, 5, 2)


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    SeqGaussSim(var₁=>param₁, var₂=>param₂, ...)

A polyalgorithm sequential Gaussian simulation solver.

Each pair `var=>param` specifies the `SeqGaussSimParam`
`param` for the simulation variable `var`. In order to avoid
boilerplate code, the constructor expects pairs of `Symbol`
and `NamedTuple` instead. See [`Kriging`](@ref) documentation
for examples.

## Parameters

* `variogram` - Variogram model (default to `GaussianVariogram()`)
* `mean`      - Simple Kriging mean
* `degree`    - Universal Kriging degree
* `drifts`    - External Drift Kriging drift functions

Latter options override former options. For example, by specifying
`drifts`, the user is telling the algorithm to ignore `degree` and
`mean`. If no option is specified, Ordinary Kriging is used by
default with the `variogram` only.

* `path`         - Simulation path (default to :random)
* `neighradius`  - Radius of search neighborhood (default to 10.)
* `maxneighbors` - Maximum number of neighbors (default to 10)
"""
@simsolver SeqGaussSim begin
  @param variogram = GaussianVariogram()
  @param mean = nothing
  @param degree = nothing
  @param drifts = nothing
  @param path = :random
  @param neighradius = 10.
  @param maxneighbors = 10
end

function preprocess(problem::SimulationProblem, solver::SeqGaussSim)
  @warn "SeqGaussSim not fully implemented, sorry!"
  # retrieve problem info
  pdata = data(problem)
  pdomain = domain(problem)

  # determine coordinate type
  T = hasdata(problem) ? coordtype(pdata) : coordtype(pdomain)

  # result of preprocessing
  preproc = Dict{Symbol,Tuple}()

  for (var, V) in variables(problem)
    # get user parameters
    if var ∈ keys(solver.params)
      varparams = solver.params[var]
    else
      varparams = SeqGaussSimParam()
    end

    # determine which Kriging variant to use
    if varparams.drifts ≠ nothing
      estimator = ExternalDriftKriging{T,V}(varaparams.variogram, varparams.drifts)
    elseif varparams.degree ≠ nothing
      estimator = UniversalKriging{T,V}(varparams.variogram, varparams.degree)
    elseif varparams.mean ≠ nothing
      estimator = SimpleKriging{T,V}(varparams.variogram, varparams.mean)
    else
      estimator = OrdinaryKriging{T,V}(varparams.variogram)
    end

    # determine which path to use
    if varparams.path == :simple
      path = SimplePath(pdomain)
    elseif varparams.path == :random
      path = RandomPath(pdomain)
    else
      error("invalid path type")
    end

    # determine which neighborhood to use
    neighborhood = BallNeighborhood(pdomain, varparams.neighradius)

    # determine maximum number of conditioning neighbors
    maxneighbors = varparams.maxneighbors

    preproc[var] = (estimator, path, neighborhood, maxneighbors)
  end

  preproc
end

function solve_single(problem::SimulationProblem, var::Symbol,
                      solver::SeqGaussSim, preproc)
  # retrieve problem info
  pdata = data(problem)
  pdomain = domain(problem)

  # unpack preprocessed parameters
  estimator, path, neighborhood, maxneighbors = preproc[var]

  # determine coordinate type
  T = hasdata(problem) ? coordtype(pdata) : coordtype(pdomain)

  # determine value type
  V = variables(problem)[var]

  # result for variable
  realization = Vector{V}(undef, npoints(pdomain))

  # keep track of simulated locations
  simulated = falses(npoints(pdomain))

  # consider data locations as already simulated
  for (loc, datloc) in datamap(problem, var)
    realization[loc] = value(pdata, datloc, var)
    simulated[loc] = true
  end

  # pre-allocate memory for coordinates
  coords = MVector{ndims(pdomain),coordtype(pdomain)}(undef)

  # simulation loop
  for location in path
    if !simulated[location]
      # find neighbors
      neighbors = neighborhood(location)

      # neighbors with previously simulated values
      neighbors = neighbors[view(simulated, neighbors)]

      # sample a subset of neighbors for computational purposes
      if length(neighbors) > maxneighbors
        neighbors = sample(neighbors, maxneighbors, replace=false)
      end

      # choose between marginal and conditional distribution
      if isempty(neighbors)
        # draw from marginal
        realization[location] = randn(V)
      else
        # build coordinates and observation arrays
        X = Matrix{T}(undef, ndims(pdomain), length(neighbors))
        for (j, neighbor) in enumerate(neighbors)
          coordinates!(X[:,j], pdomain, neighbor)
        end
        z = view(realization, neighbors)

        # build Kriging system
        fit!(estimator, X, z)

        try
          # estimate mean and variance
          coordinates!(coords, pdomain, location)
          μ, σ² = estimate(estimator, coords)

          # TODO: this portion of the code will be rewritten after
          # Julia v0.7. At the moment the linear algebra components
          # of the language do not provide a consistent, exception-free
          # way of checking for failure

          # fix possible numerical issues
          O = zero(typeof(σ²))
          σ² < O && (σ² = O)

          # draw from conditional
          realization[location] = μ + √σ²*randn(V)
        catch e
          if e isa SingularException
            # draw from marginal
            realization[location] = randn(V)
          else
            rethrow(e)
          end
        end
      end

      # mark location as simulated and continue
      simulated[location] = true
    end
  end

  realization
end


text:
"""
    random_unit_vector(N)

Generates N random unit vectors, i.e. random points on a unit sphere.
"""
function random_unit_vector(N::Int64)
    phis = 2 * pi * rand(Float64, N)
    cts = 2 * rand(Float64, N) .- 1
    sts = sqrt.(1. .- cts .* cts) # sin(acos(cts)) # max(0., )

    [SVector{3, Float64}(
        sts[i] .* cos.(phis[i]),
        sts[i] .* sin.(phis[i]),
        cts[i]
    ) for i in 1:N]
end


text:
export unsafe_swap!, rot, plshift, prshift, arshift
function unsafe_swap!(tape, i, j)
    @inbounds temp = tape[i]
    @inbounds tape[i] = tape[j]
    @inbounds tape[j] = temp
end

function rot(a, b, θ)
    c = cos(θ)
    s = sin(θ)
    a*c-b*s, a*s+b*c
end
plshift(x, n) = (x << n) | (x >> (sizeof(x)*8-n))
prshift(x, n) = (x >> n) | (x << (sizeof(x)*8-n))
"""right shift, sign extending."""
arshift(x::T, n) where T = (x >> n) | (x & (T(1) << (sizeof(x)*8-1)))


text:
using JuMP, GLPK, PyPlot, Printf , vOptGeneric
include("crea_read.jl") # permet de créer les fichiers : données bi-obj, solution Y_N, solutions relachées
include("xp_num.jl")
# https://github.com/fuofuo-gg/TER

#const GUROBI_ENV = Gurobi.Env()

# ==============================================================================
#= 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)
=#

mutable struct point
    x::Float64
    y::Float64
end

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

    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

# ----------------------------------------------------------------
# Model solutions optimales SPA bi obj
function SPA_vopt(nbvar::Int, nbcontraintes::Int, L::Array{Int,2}, c1::Array{Int,1}, c2::Array{Int,1})
    model = vModel(with_optimizer(GLPK.Optimizer))
    JuMP.set_silent(model)

    @variable(model, x[1:nbvar], Bin)

    @constraint(model, [i=1:nbcontraintes],(sum((x[j]*L[i,j]) for j in 1:nbvar)) == 1)

    @addobjective(model, Min, sum(c1[i]*x[i] for i in 1:nbvar))
    @addobjective(model, Min, sum(c2[i]*x[i] for i in 1:nbvar))

    vSolve(model, method=:epsilon, step = 0.5)

    Y_N = getY_N(model)

    return Y_N
end

# ==============================================================================
# Model solution relaxation linéaire sur un objectif
function relaxLinXG(nbvar::Int, nbcontraintes::Int, L::Array{Int,2}, c1::Array{Int,1}, c2::Array{Int,1}, vobj, delta, obj)

    model = Model(with_optimizer(GLPK.Optimizer))
#    model = Model(with_optimizer(Gurobi.Optimizer, GUROBI_ENV))
#    set_optimizer_attributes(model, "OutputFlag" => 0)

    @variable(model, 0.0 <= x[1:nbvar] <= 1.0 )
    @constraint(model, [i=1:nbcontraintes],(sum((x[j]*L[i,j]) for j in 1:nbvar)) == 1)
    if obj == 1
        @objective(model, Min, sum((c1[i])*x[i] for i in 1:nbvar))
        @constraint(model, sum((c2[i])*x[i] for i in 1:nbvar) <= vobj-delta)
    else
        @objective(model, Min, sum((c2[i])*x[i] for i in 1:nbvar))
        @constraint(model, sum((c1[i])*x[i] for i in 1:nbvar) <= vobj-delta)
    end
    optimize!(model)
    return objective_value(model), value.(x)
end

# ==============================================================================
# Model solution relaxation linéaire sur un objectif vertical
function relaxLinXG4(nbvar::Int, nbcontraintes::Int, L::Array{Int,2}, c1::Array{Int,1}, c2::Array{Int,1}, epsilon, obj)

    model = Model(with_optimizer(GLPK.Optimizer))
#    model = Model(with_optimizer(Gurobi.Optimizer, GUROBI_ENV))
#    set_optimizer_attributes(model, "OutputFlag" => 0)

    @variable(model, 0.0 <= x[1:nbvar] <= 1.0 )
    @constraint(model, [i=1:nbcontraintes],(sum((x[j]*L[i,j]) for j in 1:nbvar)) == 1)
    if obj == 1
        @objective(model, Min, sum((c1[i])*x[i] for i in 1:nbvar))
        @constraint(model, sum((c2[i])*x[i] for i in 1:nbvar) <= epsilon)
    else
        @objective(model, Min, sum((c2[i])*x[i] for i in 1:nbvar))
        @constraint(model, sum((c1[i])*x[i] for i in 1:nbvar) <= epsilon)
    end
    optimize!(model)
    return objective_value(model), value.(x)
end

# ==============================================================================
# Elabore 2 ensembles d'indices selon que xTilde[i] vaut 0 ou 1
function splitXG(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

#xTilde = [0,0,1,1,1,0,1]

# ==============================================================================
# Projete xTilde sur le polyedre X
function Δ2(L::Array{Int,2}, xTilde::Array{Int,1})

    nbcontraintes = size(L,1)
    nbvar = size(L,2)
    idxTilde0, idxTilde1 = splitXG(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:nbcontraintes],(sum((x[j]*L[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 x[i] > 0.001 && x[i] < 0.999
            return false
        end
        i+=1
    end
    return admissible
end

# ==============================================================================
# calcule la performance d'une solution 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

# ==============================================================================
# The gravity machine (Man of Steel) -> to terraform the world
#   e-contrainte avec z1 comme fonction a minimiser
#   arrondi qui evite de conduire vers un point domine par le point issu de la RL

# structure corresponding to one point generator
mutable struct tGenerateur
    xRlx :: Vector{Float64} # relaxed solution x
    yRlx :: Vector{Float64} # point corresponding to the relaxed solution x
    #
    xInt :: Vector{Int64}   # integer solution x
    yInt :: Vector{Int64}   # point corresponding to the integer solution x
    #
    xPrj :: Vector{Float64} # projected solution x
    yPrj :: Vector{Float64} # point corresponding to the projected solution x
    #
    xFea :: Bool            # indicate if a solution x is feasible or not
end

mutable struct tEnsBornant
    xEBD :: Vector{Float64} # solution x member of the dual bound set
    yEBD :: Vector{Float64} # performance y member of the dual bound set
end

mutable struct tPoint
    x :: Vector{Float64} # vector solution x (1..n)
    y :: Vector{Float64} # vector performances y (1..p)
end

# ============================================================================
#= On effectue une vérification si tous les fichiers de donnees ont bien ete
 crees avant d'executer gravity machine
=#

function verification(s::String)
    if !isfile("SPA_Databio/bio"*s*".txt")
        println("Création du fichier bi-objectif : ")
        dirname(pwd())
        SPA_creation("$s")
        println("Fichier bi-objectif cree : ")
    end
    nbvar, nbcontraintes, L, c1, c2 = SPA_loaddemandebi("SPA_Databio/bio"*s*".txt")
    if !isfile("Archive/Y_N_"*s*".txt")
        println("Création du fichier Y_N_ : ")
        start = time()
        Y_N = SPA_vopt(nbvar, nbcontraintes, L, c1, c2) #solver vopt
        temps = time() - start
        archiver_Y_N(s, Y_N, temps)
        println("Fichier Y_N_ créé : ")
    end

end

# ----------------------------------------------------------------
# Nettoyage de la valeur de vg[j].xRlx et calcul du point bi-objectif
function calcNoInteger(vg::Array{tGenerateur,1}, nbvar::Int, j::Int)
    for i in 1:nbvar
        if round(vg[j].xRlx[i], digits=3) == 0.0
            vg[j].xRlx[i] = 0.0
            vg[j].xInt[i] = 0
        elseif round(vg[j].xRlx[i], digits=3) == 1.0
            vg[j].xRlx[i] = 1.0
            vg[j].xInt[i] = 1
        else
            vg[j].xRlx[i] = round(vg[j].xRlx[i], digits=3)
            vg[j].xInt[i] = -1 # valeur sentinelle quand non entier
        end
    end
    return vg
end

# ----------------------------------------------------------------
# Selectionne les points pour constituer le cone d'interet centre sur le point generateur O
function selectPoints(j::Int, tailleSampling::Int, eb::Array{tEnsBornant,1})
    O=point(eb[j].yEBD[1], eb[j].yEBD[2])

    #choix de B
    flag = j
    while flag>1 && (eb[flag-1].yEBD[1] > eb[j].yEBD[1])
        flag -= 1
    end
    if flag == 1
        B=point(eb[1].yEBD[1], eb[1].yEBD[2]+1.0)
    else
        B=point(eb[flag-1].yEBD[1], eb[flag-1].yEBD[2])
    end

    #choix de A
    flag = j
    while flag<tailleSampling && (eb[flag+1].yEBD[1] < eb[j].yEBD[1])
        flag += 1
    end
    if flag == tailleSampling
        A=point(eb[end].yEBD[1]+1.0, eb[end].yEBD[2])
    else
        A=point(eb[flag+1].yEBD[1], eb[flag+1].yEBD[2])
    end
    return B, O, A
end

# ----------------------------------------------------------------
# Operations pour arrondir les valeurs non-entieres d'un générateur
function roundOneCone(j::Int, nbvar::Int, vg::Array{tGenerateur,1}, c1::Array{Int64,1}, c2::Array{Int64,1}, O::point, A::point, B::point)
    vg[j].yInt[1] = 0
    vg[j].yInt[2] = 0
    z1 = vg[j].yRlx[1]
    z2 = vg[j].yRlx[2]

    nbVarNonEntiere=0

    for i in 1:nbvar
        if vg[j].xInt[i] == -1
            # la variable i est non entiere
            nbVarNonEntiere += 1
            M = point( z1 - (vg[j].xRlx[i] * c1[i]) , z2 - (vg[j].xRlx[i] * c2[i]))
            #print(M)
            if !InSector(M, O, A, B)
                # le point M obtenu est hors du cone => variable i fixee a 1
                vg[j].xInt[i] = 1
                z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
            else
                # le point M obtenu est dans le cone => variable i fixee a 0
                vg[j].xInt[i] = 0
                z1 = z1 - (vg[j].xRlx[i] * c1[i])
                z2 = z2 - (vg[j].xRlx[i] * c2[i])
            end
        end
        vg[j].yInt[1] += vg[j].xInt[i] * c1[i]
        vg[j].yInt[2] += vg[j].xInt[i] * c2[i]
    end

    return vg, nbVarNonEntiere
end

# ----------------------------------------------------------------
# Operations pour arrondir les valeurs non-entieres de la projection
function roundTwoCones(j::Int, nbvar::Int, vg::Array{tGenerateur,1}, c1::Array{Int64,1}, c2::Array{Int64,1}, O::point, A::point, B::point, OO::point, param_round2::Int)
    vg[j].yInt[1] = 0
    vg[j].yInt[2] = 0
    Xapx=(Float64)[]; Yapx=(Float64)[]
    z1 = vg[j].yRlx[1]
    z2 = vg[j].yRlx[2]
    nbVarNonEntiere=0

    if param_round2 == 1
        for i in 1:nbvar
            if vg[j].xInt[i] == -1
                M = point( z1 - vg[j].xRlx[i] * c1[i] , z2 - vg[j].xRlx[i] * c2[i])
                push!(Xapx,M.x)
                push!(Yapx,M.y)
                # la variable i est non entiere
                nbVarNonEntiere += 1
                if !InSector(M, O, A, B)
                    # le point M obtenu est hors du cone => variable i fixee a 1
                    vg[j].xInt[i] = 1
                    z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                    z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                else
                    if !InSector(M, OO, B, A)
                        MM = point( z1 - vg[j].xRlx[i] * c1[i] + c1[i] , z2 - vg[j].xRlx[i] * c2[i] + c2[i])
                        if !InSector(MM, OO, B, A)
                            # le point MM obtenu est hors du cone => choix de 0 ou 1 avec la plus courte distance à O et OO
                            distm = sqrt((O.x-M.x)^2+(O.y-M.y)^2)+sqrt((OO.x-M.x)^2+(OO.y-M.y)^2)
                            dismm = sqrt((O.x-MM.x)^2+(O.y-MM.y)^2)+sqrt((OO.x-MM.x)^2+(OO.y-MM.y)^2)
                            if distm < dismm
                                vg[j].xInt[i] = 0
                                z1 = z1 - (vg[j].xRlx[i] * c1[i])
                                z2 = z2 - (vg[j].xRlx[i] * c2[i])
                            else
                                vg[j].xInt[i] = 1
                                z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                                z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                            end
                        else
                            vg[j].xInt[i] = 1
                            z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                            z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                        end
                        MM = point( z1 - (vg[j].xRlx[i] * c1[i]) + c1[i] , z2 - (vg[j].xRlx[i] * c2[i]) + c2[i])
                    else
                        # le point M obtenu est dans du cone => variable i fixee a 0
                        vg[j].xInt[i] = 0
                        #print(" ",i)
                        z1 = z1 - (vg[j].xRlx[i] * c1[i])
                        z2 = z2 - (vg[j].xRlx[i] * c2[i])
                    end
                end
            end
            vg[j].yInt[1] += vg[j].xInt[i] * c1[i]
            vg[j].yInt[2] += vg[j].xInt[i] * c2[i]
        end
    elseif param_round2 == 2
        i = 1
        while i <= nbvar
            if vg[j].xInt[i] == -1
                nbVarNonEntiere += 1
                k = i + 1
                while k <= nbvar
                    if vg[j].xInt[k] == -1
                        nbVarNonEntiere += 1
                        MOO = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] , z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k])
                        MOI = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] + c1[k] , z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k] + c2[k])
                        MIO = point( z1 - vg[j].xRlx[i] * c1[i] + c1[i] - vg[j].xRlx[k] * c1[k] , z2 - vg[j].xRlx[i] * c2[i] + c2[i] - vg[j].xRlx[k] * c2[k])
                        MII = point( z1 - vg[j].xRlx[i] * c1[i] + c1[i] - vg[j].xRlx[k] * c1[k] + c1[k] , z2 - vg[j].xRlx[i] * c2[i] + c2[i] - vg[j].xRlx[k] * c2[k] + c2[k])
                        if InSector(MOO, O, A, B) && InSector(MOO, OO, B, A)
                            vg[j].xInt[i] = 0
                            z1 = z1 - (vg[j].xRlx[i] * c1[i])
                            z2 = z2 - (vg[j].xRlx[i] * c2[i])
                            vg[j].xInt[k] = 0
                            z1 = z1 - (vg[j].xRlx[k] * c1[k])
                            z2 = z2 - (vg[j].xRlx[k] * c2[k])
                        elseif InSector(MOI, O, A, B) && InSector(MOI, OO, B, A)
                            vg[j].xInt[i] = 0
                            z1 = z1 - (vg[j].xRlx[i] * c1[i])
                            z2 = z2 - (vg[j].xRlx[i] * c2[i])
                            vg[j].xInt[k] = 1
                            z1 = z1 - (vg[j].xRlx[k] * c1[k]) + c1[k]
                            z2 = z2 - (vg[j].xRlx[k] * c2[k]) + c2[k]
                        elseif InSector(MIO, O, A, B) && InSector(MIO, OO, B, A)
                            vg[j].xInt[i] = 1
                            z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                            z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                            vg[j].xInt[k] = 0
                            z1 = z1 - (vg[j].xRlx[k] * c1[k])
                            z2 = z2 - (vg[j].xRlx[k] * c2[k])
                        elseif InSector(MII, O, A, B) && InSector(MII, OO, B, A)
                            vg[j].xInt[i] = 1
                            z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                            z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                            vg[j].xInt[k] = 1
                            z1 = z1 - (vg[j].xRlx[k] * c1[k]) + c1[k]
                            z2 = z2 - (vg[j].xRlx[k] * c2[k]) + c2[k]
                        else
                            distMOO = sqrt((O.x-MOO.x)^2+(O.y-MOO.y)^2)+sqrt((OO.x-MOO.x)^2+(OO.y-MOO.y)^2)
                            distMOI = sqrt((O.x-MOI.x)^2+(O.y-MOI.y)^2)+sqrt((OO.x-MOI.x)^2+(OO.y-MOI.y)^2)
                            distMIO = sqrt((O.x-MIO.x)^2+(O.y-MIO.y)^2)+sqrt((OO.x-MIO.x)^2+(OO.y-MIO.y)^2)
                            distMII = sqrt((O.x-MII.x)^2+(O.y-MII.y)^2)+sqrt((OO.x-MII.x)^2+(OO.y-MII.y)^2)
                            if InSector(MOO, O, A, B) && distMOO <= distMOI && distMOO <= distMIO && distMOO <= distMII
                                vg[j].xInt[i] = 0
                                z1 = z1 - (vg[j].xRlx[i] * c1[i])
                                z2 = z2 - (vg[j].xRlx[i] * c2[i])
                                vg[j].xInt[k] = 0
                                z1 = z1 - (vg[j].xRlx[k] * c1[k])
                                z2 = z2 - (vg[j].xRlx[k] * c2[k])
                            elseif InSector(MOI, O, A, B) && distMOI <= distMIO && distMOI <= distMII
                                vg[j].xInt[i] = 0
                                z1 = z1 - (vg[j].xRlx[i] * c1[i])
                                z2 = z2 - (vg[j].xRlx[i] * c2[i])
                                vg[j].xInt[k] = 1
                                z1 = z1 - (vg[j].xRlx[k] * c1[k]) + c1[k]
                                z2 = z2 - (vg[j].xRlx[k] * c2[k]) + c2[k]
                            elseif InSector(MIO, O, A, B) && distMIO <= distMII
                                vg[j].xInt[i] = 1
                                z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                                z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                                vg[j].xInt[k] = 0
                                z1 = z1 - (vg[j].xRlx[k] * c1[k])
                                z2 = z2 - (vg[j].xRlx[k] * c2[k])
                            else
                                vg[j].xInt[i] = 1
                                z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                                z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                                vg[j].xInt[k] = 1
                                z1 = z1 - (vg[j].xRlx[k] * c1[k]) + c1[k]
                                z2 = z2 - (vg[j].xRlx[k] * c2[k]) + c2[k]
                            end
                        end
                        k += nbvar + 1
                    end
                    k += 1
                end
                if k == nbvar + 1
                    nbVarNonEntiere += 1
                    MO = point( z1 - vg[j].xRlx[i] * c1[i] , z2 - vg[j].xRlx[i] * c2[i])
                    MI = point( z1 - vg[j].xRlx[i] * c1[i] + c1[i] , z2 - vg[j].xRlx[i] * c2[i] + c2[i])
                    if InSector(MO, O, A, B) && InSector(MO, OO, B, A)
                        vg[j].xInt[i] = 0
                        z1 = z1 - (vg[j].xRlx[i] * c1[i])
                        z2 = z2 - (vg[j].xRlx[i] * c2[i])
                    elseif InSector(MI, O, A, B) && InSector(MI, OO, B, A)
                        vg[j].xInt[i] = 1
                        z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                        z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                    else
                        distMO = sqrt((O.x-MO.x)^2+(O.y-MO.y)^2)+sqrt((OO.x-MO.x)^2+(OO.y-MO.y)^2)
                        distMI = sqrt((O.x-MI.x)^2+(O.y-MI.y)^2)+sqrt((OO.x-MI.x)^2+(OO.y-MI.y)^2)
                        if InSector(MO, O, A, B) && distMO < distMI
                            vg[j].xInt[i] = 0
                            z1 = z1 - (vg[j].xRlx[i] * c1[i])
                            z2 = z2 - (vg[j].xRlx[i] * c2[i])
                        else
                            vg[j].xInt[i] = 1
                            z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                            z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                        end
                    end
                end
            end
            vg[j].yInt[1] += vg[j].xInt[i] * c1[i]
            vg[j].yInt[2] += vg[j].xInt[i] * c2[i]
            i += 1
        end
    elseif param_round2 == 3
        for i in 1:nbvar
            if vg[j].xInt[i] == -1
                M = point( z1 - vg[j].xRlx[i] * c1[i] , z2 - vg[j].xRlx[i] * c2[i])
                MM = point( z1 - vg[j].xRlx[i] * c1[i] + c1[i] , z2 - vg[j].xRlx[i] * c2[i] + c2[i])
                if !InSector(M, O, A, B)
                    # le point M obtenu est hors du cone => variable i fixee a 1
                    vg[j].xInt[i] = 1
                    z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                    z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                else
                    distm = sqrt((O.x-M.x)^2+(O.y-M.y)^2)+sqrt((OO.x-M.x)^2+(OO.y-M.y)^2)
                    distmm = sqrt((O.x-MM.x)^2+(O.y-MM.y)^2)+sqrt((OO.x-MM.x)^2+(OO.y-MM.y)^2)
                    if distm < distmm
                        vg[j].xInt[i] = 0
                        z1 = z1 - (vg[j].xRlx[i] * c1[i])
                        z2 = z2 - (vg[j].xRlx[i] * c2[i])
                    else
                        vg[j].xInt[i] = 1
                        z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                        z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                    end
                end
            end
            vg[j].yInt[1] += vg[j].xInt[i] * c1[i]
            vg[j].yInt[2] += vg[j].xInt[i] * c2[i]
        end
    elseif param_round2 == 4
        i = 1
        while i <= nbvar
            if vg[j].xInt[i] == -1
                nbVarNonEntiere += 1
                k = i + 1
                while k <= nbvar
                    if vg[j].xInt[k] == -1
                        nbVarNonEntiere += 1
                        MOO = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] , z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k])
                        MOI = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] + c1[k] , z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k] + c2[k])
                        MIO = point( z1 - vg[j].xRlx[i] * c1[i] + c1[i] - vg[j].xRlx[k] * c1[k] , z2 - vg[j].xRlx[i] * c2[i] + c2[i] - vg[j].xRlx[k] * c2[k])
                        MII = point( z1 - vg[j].xRlx[i] * c1[i] + c1[i] - vg[j].xRlx[k] * c1[k] + c1[k] , z2 - vg[j].xRlx[i] * c2[i] + c2[i] - vg[j].xRlx[k] * c2[k] + c2[k])
                        distMOO = sqrt((O.x-MOO.x)^2+(O.y-MOO.y)^2)+sqrt((OO.x-MOO.x)^2+(OO.y-MOO.y)^2)
                        distMOI = sqrt((O.x-MOI.x)^2+(O.y-MOI.y)^2)+sqrt((OO.x-MOI.x)^2+(OO.y-MOI.y)^2)
                        distMIO = sqrt((O.x-MIO.x)^2+(O.y-MIO.y)^2)+sqrt((OO.x-MIO.x)^2+(OO.y-MIO.y)^2)
                        distMII = sqrt((O.x-MII.x)^2+(O.y-MII.y)^2)+sqrt((OO.x-MII.x)^2+(OO.y-MII.y)^2)
                        if InSector(MOO, O, A, B) && distMOO <= distMOI && distMOO <= distMIO && distMOO <= distMII
                            vg[j].xInt[i] = 0
                            z1 = z1 - (vg[j].xRlx[i] * c1[i])
                            z2 = z2 - (vg[j].xRlx[i] * c2[i])
                            vg[j].xInt[k] = 0
                            z1 = z1 - (vg[j].xRlx[k] * c1[k])
                            z2 = z2 - (vg[j].xRlx[k] * c2[k])
                        elseif InSector(MOI, O, A, B) && distMOI <= distMIO && distMOI <= distMII
                            vg[j].xInt[i] = 0
                            z1 = z1 - (vg[j].xRlx[i] * c1[i])
                            z2 = z2 - (vg[j].xRlx[i] * c2[i])
                            vg[j].xInt[k] = 1
                            z1 = z1 - (vg[j].xRlx[k] * c1[k]) + c1[k]
                            z2 = z2 - (vg[j].xRlx[k] * c2[k]) + c2[k]
                        elseif InSector(MIO, O, A, B) && distMIO <= distMII
                            vg[j].xInt[i] = 1
                            z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                            z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                            vg[j].xInt[k] = 0
                            z1 = z1 - (vg[j].xRlx[k] * c1[k])
                            z2 = z2 - (vg[j].xRlx[k] * c2[k])
                        else
                            vg[j].xInt[i] = 1
                            z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                            z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                            vg[j].xInt[k] = 1
                            z1 = z1 - (vg[j].xRlx[k] * c1[k]) + c1[k]
                            z2 = z2 - (vg[j].xRlx[k] * c2[k]) + c2[k]
                        end
                        k = nbvar + 1
                    end
                    k += 1
                end
                if k == nbvar + 1
                    nbVarNonEntiere += 1
                    M = point( z1 - vg[j].xRlx[i] * c1[i] , z2 - vg[j].xRlx[i] * c2[i])
                    MM = point( z1 - vg[j].xRlx[i] * c1[i] + c1[i] , z2 - vg[j].xRlx[i] * c2[i] + c2[i])
                    distm = sqrt((O.x-M.x)^2+(O.y-M.y)^2)+sqrt((OO.x-M.x)^2+(OO.y-M.y)^2)
                    distmm = sqrt((O.x-MM.x)^2+(O.y-MM.y)^2)+sqrt((OO.x-MM.x)^2+(OO.y-MM.y)^2)
                    if InSector(M, O, A, B) && distm < distmm
                        vg[j].xInt[i] = 0
                        z1 = z1 - (vg[j].xRlx[i] * c1[i])
                        z2 = z2 - (vg[j].xRlx[i] * c2[i])
                    else
                        vg[j].xInt[i] = 1
                        z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                        z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                    end
                end
            end
            vg[j].yInt[1] += vg[j].xInt[i] * c1[i]
            vg[j].yInt[2] += vg[j].xInt[i] * c2[i]
            i += 1
        end
    elseif param_round2 == 5
        i = 1
        while i <= nbvar
            if vg[j].xInt[i] == -1
                nbVarNonEntiere += 1
                k = i + 1
                while k <= nbvar
                    if vg[j].xInt[k] == -1
                        nbVarNonEntiere += 1
                        l = k + 1
                        while l <= nbvar
                            if vg[j].xInt[l] == -1
                                nbVarNonEntiere += 1
                                MOOO = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] - vg[j].xRlx[l] * c1[l],
                                              z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k] - vg[j].xRlx[l] * c2[l])
                                MOOI = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] - vg[j].xRlx[l] * c1[l] + c1[l],
                                              z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k] - vg[j].xRlx[l] * c2[l] + c2[l])
                                MOIO = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] - vg[j].xRlx[l] * c1[l] + c1[k],
                                              z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k] - vg[j].xRlx[l] * c2[l] + c2[k])
                                MOII = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] - vg[j].xRlx[l] * c1[l] + c1[k] + c1[l],
                                              z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k] - vg[j].xRlx[l] * c2[l] + c2[k] + c2[l])
                                MIOO = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] - vg[j].xRlx[l] * c1[l] + c1[i],
                                              z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k] - vg[j].xRlx[l] * c2[l] + c2[i])
                                MIOI = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] - vg[j].xRlx[l] * c1[l] + c1[i] + c1[l],
                                              z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k] - vg[j].xRlx[l] * c2[l] + c2[i] + c2[l])
                                MIIO = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] - vg[j].xRlx[l] * c1[l] + c1[i] + c1[k],
                                              z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k] - vg[j].xRlx[l] * c2[l] + c2[i] + c2[k])
                                MIII = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] - vg[j].xRlx[l] * c1[l] + c1[i] + c1[k] + c1[l],
                                              z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k] - vg[j].xRlx[l] * c2[l] + c2[i] + c2[k] + c2[l])
                                distMOOO = sqrt((O.x-MOOO.x)^2+(O.y-MOOO.y)^2)+sqrt((OO.x-MOOO.x)^2+(OO.y-MOOO.y)^2)
                                distMOOI = sqrt((O.x-MOOI.x)^2+(O.y-MOOI.y)^2)+sqrt((OO.x-MOOI.x)^2+(OO.y-MOOI.y)^2)
                                distMOIO = sqrt((O.x-MOIO.x)^2+(O.y-MOIO.y)^2)+sqrt((OO.x-MOIO.x)^2+(OO.y-MOIO.y)^2)
                                distMOII = sqrt((O.x-MOII.x)^2+(O.y-MOII.y)^2)+sqrt((OO.x-MOII.x)^2+(OO.y-MOII.y)^2)
                                distMIOO = sqrt((O.x-MIOO.x)^2+(O.y-MIOO.y)^2)+sqrt((OO.x-MIOO.x)^2+(OO.y-MIOO.y)^2)
                                distMIOI = sqrt((O.x-MIOI.x)^2+(O.y-MIOI.y)^2)+sqrt((OO.x-MIOI.x)^2+(OO.y-MIOI.y)^2)
                                distMIIO = sqrt((O.x-MIIO.x)^2+(O.y-MIIO.y)^2)+sqrt((OO.x-MIIO.x)^2+(OO.y-MIIO.y)^2)
                                distMIII = sqrt((O.x-MIII.x)^2+(O.y-MIII.y)^2)+sqrt((OO.x-MIII.x)^2+(OO.y-MIII.y)^2)
                                if InSector(MOOO, O, A, B) && distMOOO <= distMOOI && distMOOO <= distMOIO && distMOOO <= distMOII && distMOOO <= distMIOO && distMOOO <= distMIOI && distMOOO <= distMIIO && distMOOO <= distMIII
                                    vg[j].xInt[i] = 0
                                    z1 = z1 - (vg[j].xRlx[i] * c1[i])
                                    z2 = z2 - (vg[j].xRlx[i] * c2[i])
                                    vg[j].xInt[k] = 0
                                    z1 = z1 - (vg[j].xRlx[k] * c1[k])
                                    z2 = z2 - (vg[j].xRlx[k] * c2[k])
                                    vg[j].xInt[l] = 0
                                    z1 = z1 - (vg[j].xRlx[l] * c1[l])
                                    z2 = z2 - (vg[j].xRlx[l] * c2[l])
                                elseif InSector(MOOI, O, A, B) && distMOOI <= distMOIO && distMOOI <= distMOII && distMOOI <= distMIOO && distMOOI <= distMIOI && distMOOI <= distMIIO && distMOOI <= distMIII
                                    vg[j].xInt[i] = 0
                                    z1 = z1 - (vg[j].xRlx[i] * c1[i])
                                    z2 = z2 - (vg[j].xRlx[i] * c2[i])
                                    vg[j].xInt[k] = 0
                                    z1 = z1 - (vg[j].xRlx[k] * c1[k])
                                    z2 = z2 - (vg[j].xRlx[k] * c2[k])
                                    vg[j].xInt[l] = 1
                                    z1 = z1 - (vg[j].xRlx[l] * c1[l]) + c1[l]
                                    z2 = z2 - (vg[j].xRlx[l] * c2[l]) + c2[l]
                                elseif InSector(MOIO, O, A, B) && distMOIO <= distMOII && distMOIO <= distMIOO && distMOIO <= distMIOI && distMOIO <= distMIIO && distMOIO <= distMIII
                                    vg[j].xInt[i] = 0
                                    z1 = z1 - (vg[j].xRlx[i] * c1[i])
                                    z2 = z2 - (vg[j].xRlx[i] * c2[i])
                                    vg[j].xInt[k] = 1
                                    z1 = z1 - (vg[j].xRlx[k] * c1[k]) + c1[k]
                                    z2 = z2 - (vg[j].xRlx[k] * c2[k]) + c2[k]
                                    vg[j].xInt[l] = 0
                                    z1 = z1 - (vg[j].xRlx[l] * c1[l])
                                    z2 = z2 - (vg[j].xRlx[l] * c2[l])
                                elseif InSector(MOII, O, A, B) && distMOII <= distMIOO && distMOII <= distMIOI && distMOII <= distMIIO && distMOII <= distMIII
                                    vg[j].xInt[i] = 0
                                    z1 = z1 - (vg[j].xRlx[i] * c1[i])
                                    z2 = z2 - (vg[j].xRlx[i] * c2[i])
                                    vg[j].xInt[k] = 1
                                    z1 = z1 - (vg[j].xRlx[k] * c1[k]) + c1[k]
                                    z2 = z2 - (vg[j].xRlx[k] * c2[k]) + c2[k]
                                    vg[j].xInt[l] = 1
                                    z1 = z1 - (vg[j].xRlx[l] * c1[l]) + c1[l]
                                    z2 = z2 - (vg[j].xRlx[l] * c2[l]) + c2[l]
                                elseif InSector(MIOO, O, A, B) && distMIOO <= distMIOI && distMIOO <= distMIIO && distMIOO <= distMIII
                                    vg[j].xInt[i] = 1
                                    z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                                    z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                                    vg[j].xInt[k] = 0
                                    z1 = z1 - (vg[j].xRlx[k] * c1[k])
                                    z2 = z2 - (vg[j].xRlx[k] * c2[k])
                                    vg[j].xInt[l] = 0
                                    z1 = z1 - (vg[j].xRlx[l] * c1[l])
                                    z2 = z2 - (vg[j].xRlx[l] * c2[l])
                                elseif InSector(MIOI, O, A, B) && distMIOI <= distMIIO && distMIOI <= distMIII
                                    vg[j].xInt[i] = 1
                                    z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                                    z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                                    vg[j].xInt[k] = 0
                                    z1 = z1 - (vg[j].xRlx[k] * c1[k])
                                    z2 = z2 - (vg[j].xRlx[k] * c2[k])
                                    vg[j].xInt[l] = 1
                                    z1 = z1 - (vg[j].xRlx[l] * c1[l]) + c1[l]
                                    z2 = z2 - (vg[j].xRlx[l] * c2[l]) + c2[l]
                                elseif InSector(MIIO, O, A, B) && distMIIO <= distMIII
                                    vg[j].xInt[i] = 1
                                    z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                                    z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                                    vg[j].xInt[k] = 1
                                    z1 = z1 - (vg[j].xRlx[k] * c1[k]) + c1[k]
                                    z2 = z2 - (vg[j].xRlx[k] * c2[k]) + c2[k]
                                    vg[j].xInt[l] = 0
                                    z1 = z1 - (vg[j].xRlx[l] * c1[l])
                                    z2 = z2 - (vg[j].xRlx[l] * c2[l])
                                else
                                    vg[j].xInt[i] = 1
                                    z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                                    z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                                    vg[j].xInt[k] = 1
                                    z1 = z1 - (vg[j].xRlx[k] * c1[k]) + c1[k]
                                    z2 = z2 - (vg[j].xRlx[k] * c2[k]) + c2[k]
                                    vg[j].xInt[l] = 1
                                    z1 = z1 - (vg[j].xRlx[l] * c1[l]) + c1[l]
                                    z2 = z2 - (vg[j].xRlx[l] * c2[l]) + c2[l]
                                end
                                l = nbvar + 1
                                k = nbvar + 1
                            end
                            l += 1
                        end
                        if l == nbvar + 1
                            MOO = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] , z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k])
                            MOI = point( z1 - vg[j].xRlx[i] * c1[i] - vg[j].xRlx[k] * c1[k] + c1[k] , z2 - vg[j].xRlx[i] * c2[i] - vg[j].xRlx[k] * c2[k] + c2[k])
                            MIO = point( z1 - vg[j].xRlx[i] * c1[i] + c1[i] - vg[j].xRlx[k] * c1[k] , z2 - vg[j].xRlx[i] * c2[i] + c2[i] - vg[j].xRlx[k] * c2[k])
                            MII = point( z1 - vg[j].xRlx[i] * c1[i] + c1[i] - vg[j].xRlx[k] * c1[k] + c1[k] , z2 - vg[j].xRlx[i] * c2[i] + c2[i] - vg[j].xRlx[k] * c2[k] + c2[k])
                            distMOO = sqrt((O.x-MOO.x)^2+(O.y-MOO.y)^2)+sqrt((OO.x-MOO.x)^2+(OO.y-MOO.y)^2)
                            distMOI = sqrt((O.x-MOI.x)^2+(O.y-MOI.y)^2)+sqrt((OO.x-MOI.x)^2+(OO.y-MOI.y)^2)
                            distMIO = sqrt((O.x-MIO.x)^2+(O.y-MIO.y)^2)+sqrt((OO.x-MIO.x)^2+(OO.y-MIO.y)^2)
                            distMII = sqrt((O.x-MII.x)^2+(O.y-MII.y)^2)+sqrt((OO.x-MII.x)^2+(OO.y-MII.y)^2)
                            if InSector(MOO, O, A, B) && distMOO <= distMOI && distMOO <= distMIO && distMOO <= distMII
                                vg[j].xInt[i] = 0
                                z1 = z1 - (vg[j].xRlx[i] * c1[i])
                                z2 = z2 - (vg[j].xRlx[i] * c2[i])
                                vg[j].xInt[k] = 0
                                z1 = z1 - (vg[j].xRlx[k] * c1[k])
                                z2 = z2 - (vg[j].xRlx[k] * c2[k])
                            elseif InSector(MOI, O, A, B) && distMOI <= distMIO && distMOI <= distMII
                                vg[j].xInt[i] = 0
                                z1 = z1 - (vg[j].xRlx[i] * c1[i])
                                z2 = z2 - (vg[j].xRlx[i] * c2[i])
                                vg[j].xInt[k] = 1
                                z1 = z1 - (vg[j].xRlx[k] * c1[k]) + c1[k]
                                z2 = z2 - (vg[j].xRlx[k] * c2[k]) + c2[k]
                            elseif InSector(MIO, O, A, B) && distMIO <= distMII
                                vg[j].xInt[i] = 1
                                z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                                z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                                vg[j].xInt[k] = 0
                                z1 = z1 - (vg[j].xRlx[k] * c1[k])
                                z2 = z2 - (vg[j].xRlx[k] * c2[k])
                            else
                                vg[j].xInt[i] = 1
                                z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                                z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                                vg[j].xInt[k] = 1
                                z1 = z1 - (vg[j].xRlx[k] * c1[k]) + c1[k]
                                z2 = z2 - (vg[j].xRlx[k] * c2[k]) + c2[k]
                            end
                            k = nbvar + 1
                        end
                    end
                    k += 1
                end
                if k == nbvar + 1
                    M = point( z1 - vg[j].xRlx[i] * c1[i] , z2 - vg[j].xRlx[i] * c2[i])
                    MM = point( z1 - vg[j].xRlx[i] * c1[i] + c1[i] , z2 - vg[j].xRlx[i] * c2[i] + c2[i])
                    distm = sqrt((O.x-M.x)^2+(O.y-M.y)^2)+sqrt((OO.x-M.x)^2+(OO.y-M.y)^2)
                    distmm = sqrt((O.x-MM.x)^2+(O.y-MM.y)^2)+sqrt((OO.x-MM.x)^2+(OO.y-MM.y)^2)
                    if InSector(M, O, A, B) && distm < distmm
                        vg[j].xInt[i] = 0
                        z1 = z1 - (vg[j].xRlx[i] * c1[i])
                        z2 = z2 - (vg[j].xRlx[i] * c2[i])
                    else
                        vg[j].xInt[i] = 1
                        z1 = z1 - (vg[j].xRlx[i] * c1[i]) + c1[i]
                        z2 = z2 - (vg[j].xRlx[i] * c2[i]) + c2[i]
                    end
                end
            end
            vg[j].yInt[1] += vg[j].xInt[i] * c1[i]
            vg[j].yInt[2] += vg[j].xInt[i] * c2[i]
            i += 1
        end
    end

    #scatter(Xapx,Yapx, color="black", marker=".")
    return vg, nbVarNonEntiere
end

# ----------------------------------------------------------------
# Cherche l'indice de la variable |x*-x~| la plus grande
function hightest(xstar::Array{Float64, 1}, xtild::Array{Int, 1}, supr::Array{Int, 1}, nbvar::Int)
    indice = 0
    max = 0
    for i in 1:nbvar
        if i in supr
            #rien
        elseif abs(xstar[i] - xtild[i]) > max
            max = abs(xstar[i] - xtild[i])
            indice = i
        end
    end
    return indice
end

# ----------------------------------------------------------------
# Sélectionne les points pour former l'ensemble bornant primal
function kung(XFeas, YFeas)
    S = []
    for i=1:length(XFeas)
        push!(S, (XFeas[i] , YFeas[i]) )
    end
    sort!(S, by = x -> x[1])
    SN=[] ; push!(SN, S[1]) ; minYFeas = S[1][2]
    for i=2:length(XFeas)
        if S[i][2] < minYFeas
            push!(SN, S[i]) ; minYFeas = S[i][2]
        end
    end
    return SN
end

function LS(nbvar::Int, nbcontraintes::Int, L::Array{Int,2}, c1::Array{Int,1}, c2::Array{Int,1}, nb_de_un::Array{Int,1}, L_LS::Array{Int,2})
    origine = []
    desti = []
    i = 1
    while i <= nbvar-1
        j = i+1
        while j <= nbvar
            if nb_de_un[i] == nb_de_un[j]
                k = 1
            else
                k = 0
                i = j-1
                j = nbvar+1
            end
            while k > 0
                if L_LS[k,i] != L_LS[k,j]
                    k = 0
                    i = j-1
                    j = nbvar+1
                else
                    k += 1
                end
                if k > 0 && L_LS[k,i] == 0
                    if c1[i] + c2[i] > c1[j] + c2[j]
                        append!(origine, i)
                        append!(desti, j)
                        i += 1
                    else
                        append!(origine, j)
                        append!(desti, i)
                        i += 1
                    end
                    k = 0
                end
            end
            j += 1
        end
        i += 1
    end
    return origine, desti
end

function main(fname::String, tailleSampling::Int64, terraform::Int64, param_round2 = 5, param_LS = 1)

    #@printf("Running the gravity machine...\n\n")
    @assert tailleSampling>=3 "Erreur : Au moins 3 sont requis"

    verification(fname)

    # chargement de l'instance numerique ---------------------------------------
    nbvar, nbcontraintes, L, c1, c2, nb_de_un, L_LS = SPA_loaddemandebi("SPA_Databio/bio"*fname*".txt") # instance numerique de SPA
    nbobj = 2
    nbgen = tailleSampling

    start = time()

    # allocation de memoire pour la structure de donnees -----------------------
    vg=Vector{tGenerateur}(undef,nbgen)
    for j=1:nbgen
        vg[j]=tGenerateur( zeros(Float64,nbvar), zeros(Float64,nbobj),
                           zeros(Int64,nbvar),   zeros(Int64,nbobj),
                           zeros(Float64,nbvar), zeros(Float64,nbobj),
                           false
                         )
    end

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

    # listes de points pour les affichages graphiques --------------------------
    XEBD=(Float64)[]; YEBD=(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

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    #@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 = relaxLinXG(nbvar, nbcontraintes, L, c1, c2, typemax(Int), 0, 1) # opt fct 1
    z1RL1, z2RL1 = evaluerSolution(xf1RL, c1, c2)

    # calcule la valeur optimale relachee de f2 seule et le point (z1,z2) correspondant
    f2RL, xf2RL = relaxLinXG(nbvar, nbcontraintes, L, c1, c2, typemax(Int), 0, 2) # opt fct 2
    z1RL2, z2RL2 = evaluerSolution(xf2RL, c1, c2)

    #@printf("  f1_min=%8.2f ↔ f1_max=%8.2f (Δ=%.2f) \n",z1RL1, z1RL2, z1RL2-z1RL1)
    #@printf("  f2_min=%8.2f ↔ f2_max=%8.2f (Δ=%.2f) \n\n",z2RL2, z2RL1, z2RL1-z2RL2)

    if param_LS == 1
        origine, desti = LS(nbvar, nbcontraintes, L, c1, c2, nb_de_un, L_LS)
    end

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    #@printf("2) calcule les generateurs par e-contrainte avec z1 la fonction a minimiser\n\n")

    pasSampleV = 1.6*(z2RL1 - z2RL2) / ((tailleSampling-1)) # pas vertical de l'echantillonage
    pasSampleH = 1.6*(z1RL2 - z1RL1) / ((tailleSampling-1)) # pas horizontal de l'echantillo

    for j = 1:tailleSampling
        if j<floor(Int, tailleSampling/2)+1
            #@printf("  %2d : ϵ = %8.2f  ", j, z2RL1 - (j-1) * pasSampleV)

            # calcul d'une solution epsilon-contrainte -----------------------------
            fRL, eb[j].xEBD = relaxLinXG4(nbvar, nbcontraintes, L, c1, c2, z2RL1 - (j-1) * pasSampleV, 1)
            #@printf("fRL = %8.2f  ",round(fRL, digits=2))

        else
            #@printf("  %2d : ϵ = %8.2f  ", j, z1RL2 - (2 * floor(Int, tailleSampling/2) - j) * pasSampleH)

            # calcul d'une solution epsilon-contrainte -----------------------------
            fRL, eb[j].xEBD = relaxLinXG4(nbvar, nbcontraintes, L, c1, c2, z1RL2 - (2 * floor(Int, tailleSampling/2) - j) * pasSampleH, 2)
            #@printf("fRL = %8.2f  ",round(fRL, digits=2))
        end

        # nettoyage de la valeur de eb[j].xEBD et calcul du point bi-objectif --
        vg = calcNoInteger(vg, nbvar, j)

        eb[j].yEBD[1], eb[j].yEBD[2] = evaluerSolution(eb[j].xEBD, c1, c2)
        #@printf("[ %8.2f , %8.2f ] ", eb[j].yEBD[1], eb[j].yEBD[2])

        # copie de l'ensemble bornant dual dans la stru de donnees iterative ---
        vg[j].xRlx = deepcopy(eb[j].xEBD)
        vg[j].yRlx[1] = eb[j].yEBD[1]
        vg[j].yRlx[2] = eb[j].yEBD[2]
        push!(XEBD,eb[j].yEBD[1])
        push!(YEBD,eb[j].yEBD[2])

        # test d'admissibilite et marquage de la solution le cas echeant -------
        if estAdmissible(eb[j].xEBD)
            for i = 1:nbvar
                vg[j].xInt[i] = convert(Int, round(eb[j].xEBD[i]))
            end
            vg[j].yInt[1] = eb[j].yEBD[1]
            vg[j].yInt[2] = eb[j].yEBD[2]
            vg[j].xFea = true
            push!(XFeas,eb[j].yEBD[1])
            push!(YFeas,eb[j].yEBD[2])

            #ajout de la recherche locale
            if param_LS == 1
                for i in 1:length(origine)
                    if eb[j].xEBD[origine[i]] == 1

                        eb[j].xEBD[origine[i]] = 0
                        eb[j].xEBD[desti[i]] = 1

                        vg[j].yInt[1] = vg[j].yInt[1] - c1[origine[i]] + c1[desti[i]]
                        vg[j].yInt[2] = vg[j].yInt[2] - c2[origine[i]] + c2[desti[i]]
                    end
                end
                push!(XFeas, vg[j].yInt[1])
                push!(YFeas, vg[j].yInt[2])
            end

            #@printf(" Admissible \n")
        else
            vg[j].xFea = false
            #@printf(" x \n")
        end

    end
    #println("")

    # ==========================================================================
    #@printf("3) terraformation en cours...\n\n")

    for j in 1:tailleSampling
        if vg[j].xFea == false
            # Nettoyage de la valeur de vg[j].xRlx et calcul du point bi-objectif
            vg = calcNoInteger(vg, nbvar, j)

            vg[j].yRlx[1],vg[j].yRlx[2] = evaluerSolution(vg[j].xRlx, c1, c2)
            #@printf("  %2d : [ %8.2f , %8.2f ] ", j, vg[j].yRlx[1], vg[j].yRlx[2])

            # ----------------------------------------------------------------
            # Selectionne les points pour constituer le cone d'interet centre sur le point generateur O
            B, O, A = selectPoints(j, tailleSampling, eb)

            # ----------------------------------------------------------------
            # Operations pour arrondir les valeurs non-entieres de la solution relachee
            vg, nbVarNonEntiere = roundOneCone(j, nbvar, vg, c1, c2, O, A, B)

            #@printf("  → #round : %4d → [ %5d , %5d ] ", nbVarNonEntiere, vg[j].yInt[1], vg[j].yInt[2])
            push!(XInt,vg[j].yInt[1])
            push!(YInt,vg[j].yInt[2])

            trial = 0
        end

        while vg[j].xFea == false && trial < terraform
            #print(" -->")
            # ----------------------------------------------------------------
            # Projete la solution entiere sur le polytope X avec norme-L1
            fPrj, vg[j].xPrj = Δ2(L, vg[j].xInt)

            # Nettoyage de la valeur de vg[j].xPrj et calcul du point bi-objectif
            vg = calcNoInteger(vg, nbvar, j)

            vg[j].yPrj[1], vg[j].yPrj[2] = evaluerSolution(vg[j].xPrj, c1, c2)

            #@printf(" [ %8.2f , %8.2f ] ", vg[j].yPrj[1], vg[j].yPrj[2])
            push!(XProj, vg[j].yPrj[1])
            push!(YProj, vg[j].yPrj[2])

            # ----------------------------------------------------------------
            # Teste si la projection est admissible
            if estAdmissible(vg[j].xPrj)
                vg[j].xFea = true
                vg[j].yInt[1] = vg[j].yPrj[1]
                vg[j].yInt[2] = vg[j].yPrj[2]
                push!(XFeas, vg[j].yInt[1])
                push!(YFeas, vg[j].yInt[2])

                #ajout de la recherche locale
                if param_LS == 1
                    for i in 1:length(origine)
                        if vg[j].xPrj[origine[i]] == 1

                            vg[j].xPrj[origine[i]] = 0
                            vg[j].xPrj[desti[i]] = 1

                            vg[j].yInt[1] = vg[j].yInt[1] - c1[origine[i]] + c1[desti[i]]
                            vg[j].yInt[2] = vg[j].yInt[2] - c2[origine[i]] + c2[desti[i]]
                        end
                    end
                    push!(XFeas, vg[j].yInt[1])
                    push!(YFeas, vg[j].yInt[2])
                end
                #@printf(" Admissible \n")
            else
                vg[j].xFea = false
                #@printf(" x \n")
            end
            if vg[j].xFea == false && trial < (terraform - 1)

                # prepare pour l'iteration suivante
                OO = point(vg[j].yPrj[1], vg[j].yPrj[2])

                # ----------------------------------------------------------------
                # Operations pour arrondir les valeurs non-entieres de la solution relachee
                vg, nbVarNonEntiere = roundTwoCones(j, nbvar, vg, c1, c2, O, A, B, OO, param_round2)

                # ----------------------------------------------------------------
                # Gestion de 1-cycle
                if vg[j].yInt[1] == XInt[end] && vg[j].yInt[2] == YInt[end] && trial > 0
                    T = 30
                    supr = Int64[]
                    p = rand(T/2:3*T/4)
                    i = 1
                    while i <= p
                        indice = hightest(vg[j].xPrj, vg[j].xInt, supr, nbvar)
                        if indice != 0
                            append!(supr, indice)
                            if vg[j].xInt[indice] == 0
                                M = point( vg[j].yInt[1] + c1[indice] , vg[j].yInt[2] + c2[indice])
                                if InSector(M, O, A, B)
                                    if InSector(M, OO, B, A)
                                        vg[j].xInt[indice] = 1
                                        vg[j].yInt[1] += c1[indice]
                                        vg[j].yInt[2] += c2[indice]
                                    end
                                end
                            elseif vg[j].xInt[indice] == 1
                                M = point( vg[j].yInt[1] - c1[indice] , vg[j].yInt[2] - c2[indice])
                                if InSector(M, O, A, B)
                                    if InSector(M, OO, B, A)
                                        vg[j].xInt[indice] = 0
                                        vg[j].yInt[1] -= c1[indice]
                                        vg[j].yInt[2] -= c2[indice]
                                    end
                                end
                            end
                        i += 1
                        else
                            i = p+1
                        end
                    end
                end
                #@printf("  → #round : %4d → [ %5d , %5d ] ", nbVarNonEntiere, vg[j].yInt[1], vg[j].yInt[2])
                push!(XInt,vg[j].yInt[1])
                push!(YInt,vg[j].yInt[2])
            end
            trial += 1
            if length(XInt) > 1 && XInt[end] == XInt[end-1] && YInt[end] == YInt[end-1]
                trial = terraform
            end
        end
        #println("")
    end

    temps = time() - start

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

    # Donne les points relaches initiaux ---------------------------------------
    #scatter(XEBD,YEBD,color="blue", marker="x", label = "images relaxations lineaires")
    #@show XEBD
    #@show YEBD

    # Donne les points entiers -------------------------------------------------
    #scatter(XInt,YInt,color="orange", marker="s", label = " images arrondis rela lin")
    #@show XInt
    #@show YInt

    # Donne les points apres projection Δ(x,x̃) ---------------------------------
    #scatter(XProj,YProj, color="red", marker="x", label = "images projections")
    #@show XProj
    #@show YProj

    # Donne les points admissibles ---------------------------------------------
    #scatter(XFeas,YFeas, color="green", marker="o", label = "images sol admissibles")
    #@show XFeas
    #@show YFeas
    XSN = []
    YSN = []
    XSNR = []
    YSNR = []
    if length(XFeas) > 0
        SN = kung(XFeas, YFeas)
        push!(XSN, SN[1][1])
        push!(YSN, 1.1 * maximum(YInt))
        for i in 1:length(SN)-1
            push!(XSN, SN[i][1])
            push!(YSN, SN[i][2])
            push!(XSN, SN[i+1][1])
            push!(YSN, SN[i][2])
            push!(XSNR, SN[i][1])
            push!(YSNR, SN[i][2])
        end
        push!(XSN, SN[end][1])
        push!(YSN, SN[end][2])
        push!(XSN, 1.1 * maximum(XInt))
        push!(YSN, SN[end][2])
        push!(XSNR, SN[end][1])
        push!(YSNR, SN[end][2])
        #plot(XSN, YSN, color="green", marker="o", label = "ensemble bornant primal")
        #scatter(XSNR, YSNR, color="green", label = "images sol admissibles non dominées", s = 150, alpha = 0.3)
    end

    #Donne les solutions optimales du problème ---------------------------------
    X_Y_N, Y_Y_N = lecture_Y_N(fname::String)
    SM = kung(X_Y_N, Y_Y_N)
    XSM = []
    YSM = []
    XSMR = []
    YSMR = []
    for i in 1:length(SM)-1
        push!(XSM, SM[i][1])
        push!(YSM, SM[i][2])
        push!(XSM, SM[i+1][1])
        push!(YSM, SM[i][2])
        push!(XSMR, SM[i][1])
        push!(YSMR, SM[i][2])
    end
    push!(XSM, SM[end][1])
    push!(YSM, SM[end][2])
    push!(XSMR, SM[end][1])
    push!(YSMR, SM[end][2])
    #plot(XSM, YSM, color="black", label = "ensemble solutions optimales")
    #scatter(XSM, YSM, color="black", marker="+", label = "images points non dominés")
    #@show X_Y_N
    #@show Y_Y_N

    print(" & ", length(XSNR) )

    print(" & ", round(temps, digits = 2))

    println(" & ", round(area(XSM, YSM, XSN, YSN)*100, digits = 2), " \\\\")

    #legend(bbox_to_anchor=[1,1], loc=0, borderaxespad=0, fontsize = "x-small", title = "Legende")
    #PyPlot.title("f-p bi-obj conique 1 round: $fname")

end

#=
fname :
    nom de l'instance sans le .txt
tailleSampling :
    nombre de générateurs
terraform :
    nombre d'itération par générateur
Param_round2 :
    1. arrondi 1 variabale par 1 varaible avec 2 cones
    2. arrondi 2 varaibles par 2 variables avec 2 cones
    3. arrondi 1 varaible par 1 varaible avec 1 cone et 1 direction
    4. arrondi 2 varaibles par 2 varaibles avec 1 cone et 1 direction
    5. arrondi 3 varaibles par 3 varaibles avec 1 cone et 1 direction
=#

#println("Entrer nom de l'instance sans le .txt : ")
#s = chomp(readline())
#print("$s")

main("sppnw19", 30, 5)

#=
println("")
println("\\hline")
print("sppaa03 ")
main("sppaa03", 30, 5)
println("\\hline")
print("sppaa05 ")
main("sppaa05", 30, 5)
println("\\hline")
print("sppnw01 ")
main("sppnw01", 30, 5)
println("\\hline")
for i in 3:9
    b = string(i)
    s = "sppnw0"*b
    print("$s ")
    main("$s", 30, 5)
    println("\\hline")
end

for i in 19:43
    b = string(i)
    s = "sppnw"*b
    print("$s ")
    main("$s", 30, 5)
    println("\\hline")
end
=#
nothing

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
module Gcode

import Base.write

export G, write, move, meander

type G
    output::IO
    echo::Bool
    header::String
    footer::String
    position::Dict
    absolute::Bool
    feedrate::Real
    flavor::String
    view::Bool
end

function G(;output::String = "",
            echo::Bool = false,
            header::String = "",
            footer::String = "",
            position::Dict = Dict{String, Float64}(),
            absolute::Bool = true,
            feedrate::Real = 0,
            flavor::String = "Marlin",
            view::Bool = false)

    io = (output != "") ? open(output, "w") : STDOUT
    if header != ""
        write(io, readall(open(header, "r")))
    end

    return G(io,echo,header,footer,position,absolute,feedrate,flavor,view)
end

function format_args(;args...)
    str = ""
    for (key, val) in args
        str *= uppercase(string(" ", key, val))
    end
    return str
end

function write(g::G, x)
    str = string(x)
    write(g.output, str*"\n")
    if g.echo; println(str); end
end

function set_home!(g::G; args...):
    """ Set the current position to the given position without moving.

Example
-------
>>> # set the current position to X=0, Y=0
>>> g.set_home(0, 0)

"""
    write(g, "G92"*format_args(;args...))

    updateposition!(g, args)
end

function relative!(g::G):
    """ Enter relative movement mode, in general this method should not be
used, most methods handle it automatically.

"""
    write(g, "G91")
    g.absolute = false
end

function absolute!(g::G):
    """ Enter absolute movement mode, in general this method should not be
used, most methods handle it automatically.

"""
    write(g,"G90")
    g.absolute = true
end

function update_position!(g::G; args...)
    for (axis, delta) in args
        axis_str = string(axis)
        if haskey(g.position, axis_str)
            if !g.absolute
                g.position[axis_str] += delta
            else
                g.position[string(axis)] = delta
            end
        else
            merge!(g.position, [axis_str=>delta])
        end
    end
end

function feedrate!(g::G, rate)
    """ Set the feed rate (tool head speed) in mm/s

Parameters
----------
rate : float
The speed to move the tool head in mm/s.

"""
    write(g, rate)
end

function dwell(g::G, time)
    """ Pause code executions for the given amount of time.

Parameters
----------
time : float
Time in seconds to pause code execution.

"""
    write(g, "G4 P"*string(time))
end

function teardown(g::G)
    """ Close the outfile file after writing the footer if opened. This
method must be called once after all commands.

"""
    if g.footer != ""
        write(g.output, readall(open(g.footer, "r")))
    end
end

function home(self):
    """ Move the tool head to the home position (X=0, Y=0).
"""
    self.abs_move(x=0, y=0)
end

function move(g::G ; args...)
    """ Move the tool head to the given position. This method operates in
relative mode unless a manual call to `absolute` was given previously.
If an absolute movement is desired, the `abs_move` method is
recommended instead.

Examples
--------
>>> # move the tool head 10 mm in x and 10 mm in y
>>> g.move(x=10, y=10)
>>> # the x and y keywords may be omitted:
>>> g.move(10, 10)

>>> # move the A axis up 20 mm
>>> g.move(A=20)

"""
    update_position!(g;args...)

    args = format_args(;args...)
    write(g, "G1 " * args)
end

function abs_move(g::G ; args...)
    """ Same as `move` method, but positions are interpreted as absolute.
"""
    if !g.absolute
        absolute!(g)
        move(g;args...)
        relative!(g)
    else #already absolute
        move(g;args...)
    end
end

#function arc(self, direction='CW', radius='auto', helix_dim=None, helix_len=0,
#        **kwargs):
#    """ Arc to the given point with the given radius and in the given
#direction. If helix_dim and helix_len are specified then the tool head
#will also perform a linear movement through the given dimension while
#completing the arc.

#Parameters
#----------
#points : floats
#Must specify two points as kwargs, e.g. X=5, Y=5
#direction : str (either 'CW' or 'CCW') (default: 'CW')
#The direction to execute the arc in.
#radius : 'auto' or float (default: 'auto')
#The radius of the arc. A negative value will select the longer of
#the two possible arc segments. If auto is selected the radius will
#be set to half the linear distance to desired point.
#helix_dim : str or None (default: None)
#The linear dimension to complete the helix through
#helix_len : float
#The length to move in the linear helix dimension.

#Examples
#--------
#>>> # arc 10 mm up in y and 10 mm over in x with a radius of 20.
#>>> g.arc(x-10, y=10, radius=20)

#>>> # move 10 mm up on the A axis, arcing through y with a radius of 20
#>>> g.arc(A=10, y=0, radius=20)

#>>> # arc through x and y while moving linearly on axis A
#>>> g.arc(x=10, y=10, radius=50, helix_dim='A', helix_len=5)

#"""
#    msg = 'Must specify point with 2 dimensions as keywords, e.g. X=0, Y=10'
#    if len(kwargs) != 2:
#        raise RuntimeError(msg)
#    dimensions = [k.lower() for k in kwargs.keys()]
#    if 'x' in dimensions and 'y' in dimensions:
#        plane_selector = 'G17' # XY plane
#        axis = helix_dim
#    elif 'x' in dimensions:
#        plane_selector = 'G18' # XZ plane
#        dimensions.remove('x')
#        axis = dimensions[0].upper()
#    elif 'y' in dimensions:
#        plane_selector = 'G19' # YZ plane
#        dimensions.remove('y')
#        axis = dimensions[0].upper()
#    else:
#        raise RuntimeError(msg)

#    if direction == 'CW':
#        command = 'G2'
#    elif direction == 'CCW':
#        command = 'G3'

#    values = kwargs.values()
#    if self.movement_mode == 'relative':
#        dist = math.sqrt(values[0] ** 2 + values[1] ** 2)
#    else:
#        k = kwargs.keys()
#        cp = self.current_position
#        dist = math.sqrt(
#            (cp[k[0]] - values[0]) ** 2 + (cp[k[1]] - values[1]) ** 2
#        )
#    if radius == 'auto':
#        radius = dist / 2.0
#    elif radius < dist / 2.0:
#        msg = 'Radius {} to small for distance {}'.format(radius, dist)
#        raise RuntimeError(msg)

#    if axis is not None:
#        self.write('G16 X Y {}'.format(axis)) # coordinate axis assignment
#    self.write(plane_selector)
#    args = ' '.join([(k.upper() + str(v)) for k, v in kwargs.items()])
#    if helix_dim is None:
#        self.write('{} {} R{}'.format(command, args, radius))
#    else:
#        self.write('{} {} R{} G1 {}{}'.format(command, args, radius,
#                                              helix_dim.upper(), helix_len))
#        kwargs[helix_dim] = helix_len

#    self._update_current_position(**kwargs)
#end

#function abs_arc(self, direction='CW', radius='auto', **kwargs):
#    """ Same as `arc` method, but positions are interpreted as absolute.
#"""
#    self.absolute()
#    self.arc(direction=direction, radius=radius, **kwargs)
#    self.relative()
#end

function rect(g, x, y; direction::String="CW", start::String="LL"):
    """ Trace a rectangle with the given width and height.

Parameters
----------
x : float
The width of the rectangle in the x dimension.
y : float
The heigh of the rectangle in the y dimension.
direction : str (either 'CW' or 'CCW') (default: 'CW')
Which direction to complete the rectangle in.
start : str (either 'LL', 'UL', 'LR', 'UR') (default: 'LL')
The start of the rectangle - L/U = lower/upper, L/R = left/right
This assumes an origin in the lower left.

Examples
--------
>>> # trace a 10x10 clockwise square, starting in the lower left corner
>>> g.rect(10, 10)

>>> # 1x5 counterclockwise rect starting in the upper right corner
>>> g.rect(1, 5, direction='CCW', start='UR')

"""
    start = uppercase(start)
    direction = uppercase(direction)
    if direction == "CW"
        if start == "LL"
            move(g, y=y)
            move(g, x=x)
            move(g, y=-y)
            move(g, x=-x)
        elseif start == "UL"
            move(g, x=x)
            move(g, y=-y)
            move(g, x=-x)
            move(g, y=y)
        elseif start == "UR"
            move(g, y=-y)
            move(g, x=-x)
            move(g, y=y)
            move(g, x=x)
        elseif start == "LR"
            move(g, x=-x)
            move(g, y=y)
            move(g, x=x)
            move(g, y=-y)
        end
    elseif direction == "CCW"
        if start == "LL"
            move(g, x=x)
            move(g, y=y)
            move(g, x=-x)
            move(g, y=-y)
        elseif start == "UL"
            move(g, y=-y)
            move(g, x=x)
            move(g, y=y)
            move(g, x=-x)
        elseif start == "UR"
            move(g, x=-x)
            move(g, y=-y)
            move(g, x=x)
            move(g, y=y)
        elseif start == "LR"
            move(g, y=y)
            move(g, x=-x)
            move(g, y=-y)
            move(g, x=x)
        end
    end
end

function meander(g::G, x, y, spacing; start="LL", orientation="x", tail=false):
    """ Infill a rectangle with a square wave meandering pattern. If the
relevant dimension is not a multiple of the spacing, the spacing will
be tweaked to ensure the dimensions work out.

Parameters
----------
x : float
The width of the rectangle in the x dimension.
y : float
The heigh of the rectangle in the y dimension.
spacing : float
The space between parallel meander lines.
start : str (either 'LL', 'UL', 'LR', 'UR') (default: 'LL')
The start of the meander - L/U = lower/upper, L/R = left/right
This assumes an origin in the lower left.
orientation : str ('x' or 'y') (default: 'x')

Examples
--------
>>> # meander through a 10x10 sqaure with a spacing of 1mm starting in
>>> # the lower left.
>>> g.meander(10, 10, 1)

>>> # 3x5 meander with a spacing of 1 and with parallel lines through y
>>> g.meander(3, 5, spacing=1, orientation='y')

>>> # 10x5 meander with a spacing of 2 starting in the upper right.
>>> g.meander(10, 5, 2, start='UR')

"""
    start = uppercase(start)
    if start == "UL"
        x, y = x, -y
    elseif start == "UR"
        x, y = -x, -y
    elseif start == "LR"
        x, y = -x, y
    end

    # Major axis is the parallel lines, minor axis is the jog.
    orientation = lowercase(orientation)
    if orientation == "x"
        major, major_name = x, :x
        minor, minor_name = y, :y
    else
        major, major_name = y, :y
        minor, minor_name = x, :x
    end

    if minor > 0
        passes = ceil(minor / spacing)
    else
        passes = abs(floor(minor / spacing))
    end
    actual_spacing = minor / passes
    if abs(actual_spacing) != spacing
        msg = "meander spacing updated from $spacing to $actual_spacing"
        warn(msg)
        write(g, ";"*msg)
    end
    spacing = actual_spacing
    sign = 1
    relative!(g)
    for _ in 1:passes
        move(g;{major_name => (sign * major)}...)
        move(g;{minor_name => spacing}...)
        sign = -1 * sign
    end
    if tail == false
        move(g;{major_name => (sign * major)}...)
    end
end

end # module


text:
# This file is part of Kpax3. License is MIT.

include("local_mode.jl")

include("selection.jl")
include("crossover.jl")
include("mutation.jl")

include("algorithm.jl")


text:
include("./produce_data.jl")
using NPZ
using HDF5

function main(;
    )
    # create data
    data = produce_data_avalanches_sandpile_illustration()
    npzwrite(@sprintf("%s/data/data_avalanche_sandpile_illustration.npz", path_project),data)
end

# call function
main()


text:
# Classic Hodgkin Huxley neuron with a "current pulse" stimulus
using Conductor, IfElse, Unitful, ModelingToolkit
import Unitful: mV, mS, cm, µm, pA, nA, mA, µA, ms
import Conductor: Na, K # shorter aliases for Sodium/Potassium

Vₘ = MembranePotential(-65mV) # set default V₀ = -65mV

nav_kinetics = [
    Gate(AlphaBeta,
         IfElse.ifelse(Vₘ == -40.0, 1.0, (0.1*(Vₘ + 40.0))/(1.0 - exp(-(Vₘ + 40.0)/10.0))),
         4.0*exp(-(Vₘ + 65.0)/18.0), 3, name = :m)
    Gate(AlphaBeta,
         0.07*exp(-(Vₘ+65.0)/20.0),
         1.0/(1.0 + exp(-(Vₘ + 35.0)/10.0)), name = :h)]

kdr_kinetics = [
    Gate(AlphaBeta,
         IfElse.ifelse(Vₘ == -55.0, 0.1, (0.01*(Vₘ + 55.0))/(1.0 - exp(-(Vₘ + 55.0)/10.0))),
         0.125 * exp(-(Vₘ + 65.0)/80.0),
         4, name = :n)]

@named NaV = IonChannel(Sodium, nav_kinetics, max_g = 120mS/cm^2)
@named Kdr = IonChannel(Potassium, kdr_kinetics, max_g = 36mS/cm^2)
@named leak = IonChannel(Leak, max_g = 0.3mS/cm^2)

channels = [NaV, Kdr, leak];
reversals = Equilibria([Na => 50.0mV, K => -77.0mV, Leak => -54.4mV])

@named Iₑ = IonCurrent(NonIonic)
electrode_pulse = Iₑ ~ IfElse.ifelse(t > 100.0, IfElse.ifelse(t < 200.0, ustrip(Float64, µA, 400pA), 0.0), 0.0)

@named neuron = CompartmentSystem(Vₘ, channels, reversals;
                                  geometry = Sphere(radius = 20µm),
                                  stimuli = [electrode_pulse])

sim = Simulation(neuron, time = 300ms)

using OrdinaryDiffEq

solution = solve(sim, Rosenbrock23())

using Plots
# Plot at 5kHz sampling
plot(solution; plotdensity=Int(1500), size=(1200,800))


text:
# FFT operators: joDFT

## helper module
module joDFT_etc
    using JOLI: jo_convert
    using FFTW
    ### planned
    function apply_fft(pln::FFTW.cFFTWPlan,v::Vector{vdt},ms::Dims,rdt::DataType,centered::Bool) where vdt<:Union{AbstractFloat,Complex}
        mp=prod(ms)
        rv=reshape(v,ms)
        rv=(pln*rv)/sqrt(mp)
        if centered rv=fftshift(rv) end
        rv=vec(rv)
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    function apply_fft(pln::FFTW.cFFTWPlan,v::Matrix{vdt},ms::Dims,rdt::DataType,centered::Bool) where vdt<:Union{AbstractFloat,Complex}
        lms=length(ms)
        nvc=size(v,2)
        msc=(ms...,nvc)
        dims=[1:lms...]
        mp=prod(ms)
        rv=reshape(v,msc)
        pf=plan_fft(rv,dims)
        rv=(pf*rv)/sqrt(mp)
        if centered rv=fftshift(rv,dims) end
        rv=reshape(rv,(mp,nvc))
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    function apply_ifft(pln::AbstractFFTs.ScaledPlan,v::Vector{vdt},ms::Dims,rdt::DataType,centered::Bool) where vdt<:Union{AbstractFloat,Complex}
        mp=prod(ms)
        rv=reshape(v,ms)
        if centered rv=ifftshift(rv) end
        rv=(pln*rv)*sqrt(mp)
        rv=vec(rv)
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    function apply_ifft(pln::AbstractFFTs.ScaledPlan,v::Matrix{vdt},ms::Dims,rdt::DataType,centered::Bool) where vdt<:Union{AbstractFloat,Complex}
        lms=length(ms)
        nvc=size(v,2)
        msc=(ms...,nvc)
        dims=[1:lms...]
        mp=prod(ms)
        rv=reshape(v,msc)
        pf=plan_ifft(rv,dims)
        if centered rv=fftshift(rv,dims) end
        rv=(pf*rv)*sqrt(mp)
        rv=reshape(rv,(prod(ms),nvc))
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    ### not planned
    function apply_fft(v::Vector{vdt},ms::Dims,rdt::DataType,centered::Bool) where vdt<:Union{AbstractFloat,Complex}
        mp=prod(ms)
        rv=reshape(v,ms)
        rv=fft(rv)/sqrt(mp)
        if centered rv=fftshift(rv) end
        rv=vec(rv)
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    function apply_fft(v::Matrix{vdt},ms::Dims,rdt::DataType,centered::Bool) where vdt<:Union{AbstractFloat,Complex}
        lms=length(ms)
        nvc=size(v,2)
        msc=(ms...,nvc)
        dims=[1:lms...]
        mp=prod(ms)
        rv=reshape(v,msc)
        pf=plan_fft(rv,dims)
        rv=(pf*rv)/sqrt(mp)
        if centered rv=fftshift(rv,dims) end
        rv=reshape(rv,(mp,nvc))
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    function apply_ifft(v::Vector{vdt},ms::Dims,rdt::DataType,centered::Bool) where vdt<:Union{AbstractFloat,Complex}
        mp=prod(ms)
        rv=reshape(v,ms)
        if centered rv=ifftshift(rv) end
        rv=ifft(rv)*sqrt(mp)
        rv=vec(rv)
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    function apply_ifft(v::Matrix{vdt},ms::Dims,rdt::DataType,centered::Bool) where vdt<:Union{AbstractFloat,Complex}
        lms=length(ms)
        nvc=size(v,2)
        msc=(ms...,nvc)
        dims=[1:lms...]
        mp=prod(ms)
        rv=reshape(v,msc)
        pf=plan_ifft(rv,dims)
        if centered rv=fftshift(rv,dims) end
        rv=(pf*rv)*sqrt(mp)
        rv=reshape(rv,(prod(ms),nvc))
        rv=jo_convert(rdt,rv,false)
        return rv
    end
end
using .joDFT_etc

# constructors
export joDFT
"""
    julia> op = joDFT(m[,n[,...]];[plannned=...,][centered=...,][DDT=...,][RDT=...,][name=...])

Multi-dimensional FFT transform over fast dimension(s)

# Signature

    function joDFT(ms::Integer...;
        planned::Bool=true, centered::Bool=false,
        DDT::DataType=joFloat, RDT::DataType=(DDT<:Real ? Complex{DDT} : DDT),
        name::String="joDFT")

# Arguments

- `m`: dimension
- optional
    - `n,...`: more dimensions
- keywords
    - `planned`: use FFTW plan
    - `centered`: produce centered coefficients
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Notes

- multidimensional image must be vectorized
- if DDT:<Real then imaginary part will be neglected for transpose/adjoint
- if you intend to use joDFT in remote* calls, you have to either set planned=false or create the operator on the worker
- joDFT is always planned if applied to multi-vector

# Examples

1D FFT

    joDFT(m)

1D FFT with centered coefficients

    joDFT(m; centered=true)

1D FFT without the precomputed plan

    joDFT(m; planned=false)

2D FFT

    joDFT(m,n)

examples with DDT/RDT

    joDFT(m; DDT=Float32)
    joDFT(m; DDT=Float32,RDT=ComplexF64)

"""
function joDFT(ms::Integer...;
    planned::Bool=true, centered::Bool=false,
    DDT::DataType=joFloat, RDT::DataType=(DDT<:Real ? Complex{DDT} : DDT),
    name::String="joDFT")

    if planned
        pf=plan_fft(zeros(ms))
        ipf=plan_ifft(zeros(ms))
        if centered
            return joLinearFunction_A(prod(ms),prod(ms),
                v1->joDFT_etc.apply_fft(pf,v1,ms,RDT,true),
                v2->joDFT_etc.apply_ifft(ipf,v2,ms,DDT,true),
                v3->joDFT_etc.apply_ifft(ipf,v3,ms,DDT,true),
                v4->joDFT_etc.apply_fft(pf,v4,ms,RDT,true),
                DDT,RDT;fMVok=true,iMVok=true,
                name=name*"_pc"
                )
        else
            return joLinearFunction_A(prod(ms),prod(ms),
                v1->joDFT_etc.apply_fft(pf,v1,ms,RDT,false),
                v2->joDFT_etc.apply_ifft(ipf,v2,ms,DDT,false),
                v3->joDFT_etc.apply_ifft(ipf,v3,ms,DDT,false),
                v4->joDFT_etc.apply_fft(pf,v4,ms,RDT,false),
                DDT,RDT;fMVok=true,iMVok=true,
                name=name*"_p"
                )
        end
    else
        if centered
            return joLinearFunction_A(prod(ms),prod(ms),
                v1->joDFT_etc.apply_fft(v1,ms,RDT,true),
                v2->joDFT_etc.apply_ifft(v2,ms,DDT,true),
                v3->joDFT_etc.apply_ifft(v3,ms,DDT,true),
                v4->joDFT_etc.apply_fft(v4,ms,RDT,true),
                DDT,RDT;fMVok=true,iMVok=true,
                name=name*"_c"
                )
        else
            return joLinearFunction_A(prod(ms),prod(ms),
                v1->joDFT_etc.apply_fft(v1,ms,RDT,false),
                v2->joDFT_etc.apply_ifft(v2,ms,DDT,false),
                v3->joDFT_etc.apply_ifft(v3,ms,DDT,false),
                v4->joDFT_etc.apply_fft(v4,ms,RDT,false),
                DDT,RDT;fMVok=true,iMVok=true,
                name=name
                )
        end
    end
end


text:
nums = [[parse(Int64, l) for l in line] for line in readlines("inputs/1-5/03.txt")]
lines = readlines("inputs/1-5/03.txt")

function day3(lines)
  nums = [[parse(UInt8, l) for l in line] for line in lines]
  nums_backup = copy(nums)
  oxygen_vec = 0
  for digit in 1:12
    Σ = sum(map(x -> x[digit], nums) / length(nums))
    signal = Σ >= 0.5 ? 1 : 0
    signal = Σ/length(nums) == .5 ? 1 : signal

    filter!(x -> x[digit]==signal,nums)

    if length(nums) == 1
      oxygen_vec = nums[1]
      break
    end
  end

  nums = nums_backup
  co2_vec = 0
  for digit in 1:12
    Σ = sum(map(x -> x[digit], nums) / length(nums))
    signal = Σ >= 0.5 ? 0 : 1
    signal = Σ/length(nums) == .5 ? 0 : signal

    filter!(x -> x[digit]==signal,nums)

    if length(nums) == 1
      co2_vec = nums[1]
      break
    end
  end
  return parse(Int64,join(oxygen_vec),base=2) * parse(Int64,join(co2_vec),base=2)
end

# println(day3(nums))
using BenchmarkTools
@btime day3(lines)

print(day3(lines))


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

struct UsageApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets the current usage count and the limit for the resources under the subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: UsageListResult
"""
function usageList(_api::UsageApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", UsageListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Storage/usages", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export usageList