text:
@enum Neighborhood KP_1_1 KP_1_2 KP_2_1 KP_0_1
@enum Dominance DOMINATES DOMINATED NOTHING

function d(solA::Solution, solB::Solution)
    return sqrt((solB.y[1] - solA.y[1])^2 + (solB.y[2] - solA.y[2])^2)
end

function d(yA::Vector{Float64}, yB::Vector{Float64})
    #@assert length(yA) == 2 && length(yB) == 2 "We only work with 2D"
    return sqrt((yB[1] - yA[1])^2 + (yB[2] - yA[2])^2)
end

function d(yA::Vector{Float64}, yB::Vector{Float64})
    #@assert length(yA) == 2 && length(yB) == 2 "We only work with 2D"
    return sqrt((yB[1] - yA[1])^2 + (yB[2] - yA[2])^2)
end

function middle(yA::Vector{Float64}, yB::Vector{Float64})
    #@assert length(yA) == 2 && length(yB) == 2 "We only work with 2D"
    return [(yA[1]+yB[1])/2, (yA[2]+yB[2])/2]
end

function evaluate(y::Vector{Float64}, consecutiveSet::Vector{PairOfSolution}; fastEvaluate = false)
    #@assert length(y) == 2 "We only work with 2D"
    minEval = Inf
    for pairOfSol in consecutiveSet
        solL,solR = pairOfSol.solL, pairOfSol.solR
        middlePoint = middle(solL.y,solR.y)
        if !fastEvaluate
            minLengthTriangle = d(solL,solR)
            nadir = [solL.y[1],solR.y[2]]
            if d(solL.y,nadir) < minLengthTriangle
                minLengthTriangle = d(solL.y,nadir)
            end
            if d(solR.y,nadir) < minLengthTriangle
                minLengthTriangle = d(solR.y,nadir)
            end
            w = 1/minLengthTriangle
            evalSol = d(y,middlePoint)*w
        else
            evalSol = d(y,middlePoint)
        end
        if evalSol < minEval
            minEval = evalSol
        end
    end
    return minEval
end

function bestKP01(sol::Solution, consecutiveSet::Vector{PairOfSolution}, prob::Problem; fastEvaluate = false)
    bestKP = [0]
    bestEval = evaluate(sol.y,consecutiveSet,fastEvaluate=fastEvaluate)
    evalBuffer = 0.0
    for up in 1:prob.nbVar
        if sol.x[up] == 0
            # test is the kp is feasible
            if sol.w + prob.weights[up] <= prob.maxWeight # feasible
                #println("feasible : $(sol.w) + $(prob.weights[up]) <= $(prob.maxWeight)")
                # evaluate new solution
                y = [sol.y[1]+prob.profits[1,up],sol.y[2]+prob.profits[2,up]]
                evalBuffer = evaluate(y,consecutiveSet,fastEvaluate=fastEvaluate)
                if evalBuffer < bestEval
                    bestEval = evalBuffer
                    bestKP[1] = up
                end
            end
        end
    end
    return bestKP
end

function firstKP01(sol::Solution, consecutiveSet::Vector{PairOfSolution}, prob::Problem; fastEvaluate = false)
    bestKP = [0]
    bestEval = evaluate(sol.y,consecutiveSet,fastEvaluate=fastEvaluate)
    evalBuffer = 0.0
    improved = false
    up = 0
    while !improved && up < prob.nbVar
        up += 1
        if sol.x[up] == 0
            # test is the kp is feasible
            if sol.w + prob.weights[up] <= prob.maxWeight # feasible
                #println("feasible : $(sol.w) + $(prob.weights[up]) <= $(prob.maxWeight)")
                # evaluate new solution
                y = [sol.y[1]+prob.profits[1,up],sol.y[2]+prob.profits[2,up]]
                evalBuffer = evaluate(y,consecutiveSet,fastEvaluate=fastEvaluate)
                if evalBuffer < bestEval
                    improved = true
                    bestEval = evalBuffer
                    bestKP[1] = up
                end
            end
        end
    end
    return bestKP
end

function bestKP12(sol::Solution, consecutiveSet::Vector{PairOfSolution}, prob::Problem; fastEvaluate = false)
    bestKP = [0,0,0]
    bestEval = evaluate(sol.y,consecutiveSet,fastEvaluate=fastEvaluate)
    evalBuffer = 0.0
    for down in 1:prob.nbVar
        if sol.x[down] == 1
            for up in 1:prob.nbVar
                if sol.x[up] == 0
                    # test if the kp could be feasible
                    if sol.w - prob.weights[down] + prob.weights[up] < prob.maxWeight
                        for up2 in 1:prob.nbVar
                            if sol.x[up2] == 0
                                # test is the kp is feasible
                                if sol.w + prob.weights[up] + prob.weights[up2] - prob.weights[down] <= prob.maxWeight # feasible
                                    # evaluate new solution
                                    y = [sol.y[1]+prob.profits[1,up]+prob.profits[1,up2]-prob.profits[1,down],sol.y[2]+prob.profits[2,up]+prob.profits[2,up2]-prob.profits[2,down]]
                                    evalBuffer = evaluate(y,consecutiveSet,fastEvaluate=fastEvaluate)
                                    if evalBuffer < bestEval
                                        bestEval = evalBuffer
                                        bestKP[1] = down
                                        bestKP[2] = up
                                        bestKP[3] = up2
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return bestKP
end

function firstKP12(sol::Solution, consecutiveSet::Vector{PairOfSolution}, prob::Problem; fastEvaluate = false)
    bestKP = [0,0,0]
    bestEval = evaluate(sol.y,consecutiveSet,fastEvaluate=fastEvaluate)
    evalBuffer = 0.0
    improved = false
    down = 0
    while !improved && down < prob.nbVar
        down += 1
        if sol.x[down] == 1
            up = 0
            while !improved && up < prob.nbVar
                up += 1
                if sol.x[up] == 0
                    # test if the kp could be feasible
                    if sol.w - prob.weights[down] + prob.weights[up] < prob.maxWeight
                        up2 = 0
                        while !improved && up2 < prob.nbVar
                            up2 += 1
                            if sol.x[up2] == 0
                                # test is the kp is feasible
                                if sol.w + prob.weights[up] + prob.weights[up2] - prob.weights[down] <= prob.maxWeight # feasible
                                    # evaluate new solution
                                    y = [sol.y[1]+prob.profits[1,up]+prob.profits[1,up2]-prob.profits[1,down],sol.y[2]+prob.profits[2,up]+prob.profits[2,up2]-prob.profits[2,down]]
                                    evalBuffer = evaluate(y,consecutiveSet,fastEvaluate=fastEvaluate)
                                    if evalBuffer < bestEval
                                        improved = true
                                        bestEval = evalBuffer
                                        bestKP[1] = down
                                        bestKP[2] = up
                                        bestKP[3] = up2
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return bestKP
end

function bestKP11(sol::Solution, consecutiveSet::Vector{PairOfSolution}, prob::Problem, tEval::Float64, tE2::Float64; fastEvaluate = false)
    startFun = time()
    bestKP = [0,0]
    bestEval = evaluate(sol.y,consecutiveSet,fastEvaluate=fastEvaluate)
    evalBuffer = 0.0
    for down in 1:prob.nbVar
        if sol.x[down] == 1
            for up in 1:prob.nbVar
                if sol.x[up] == 0
                    # test is the kp is feasible
                    if sol.w + prob.weights[up] - prob.weights[down] <= prob.maxWeight # feasible
                        # evaluate new solution
                        y = [sol.y[1]+prob.profits[1,up]-prob.profits[1,down],sol.y[2]+prob.profits[2,up]-prob.profits[2,down]]
                        start = time()
                        evalBuffer = evaluate(y,consecutiveSet,fastEvaluate=fastEvaluate)
                        tEval += time() - start
                        if evalBuffer < bestEval
                            bestEval = evalBuffer
                            bestKP[1] = down
                            bestKP[2] = up
                        end
                    end
                end
            end
        end
    end
    return bestKP, tEval, tE2 + time() - startFun
end

function firstKP11(sol::Solution, consecutiveSet::Vector{PairOfSolution}, prob::Problem; fastEvaluate = false)
    bestKP = [0,0]
    bestEval = evaluate(sol.y,consecutiveSet,fastEvaluate=fastEvaluate)
    evalBuffer = 0.0
    improved = false
    down = 0
    while !improved && down < prob.nbVar
        down += 1
        if sol.x[down] == 1
            up = 0
            while !improved && up < prob.nbVar
                up += 1
                if sol.x[up] == 0
                    if sol.w + prob.weights[up] - prob.weights[down] <= prob.maxWeight # feasible
                        # evaluate new solution
                        y = [sol.y[1]+prob.profits[1,up]-prob.profits[1,down],sol.y[2]+prob.profits[2,up]-prob.profits[2,down]]
                        evalBuffer = evaluate(y,consecutiveSet,fastEvaluate=fastEvaluate)
                        if evalBuffer < bestEval
                            improved = true
                            bestEval = evalBuffer
                            bestKP[1] = down
                            bestKP[2] = up
                        end
                    end
                end
            end
        end
    end
    return bestKP
end

function updateSolByKP(sol::Solution,kp::Vector{Int}, neighborhood::Neighborhood, prob::Problem)
    newX = copy(sol.x)
    if neighborhood == KP_1_1
        newX[kp[1]] = 0
        newX[kp[2]] = 1
        newY = [sol.y[1]+prob.profits[1,kp[2]]-prob.profits[1,kp[1]],sol.y[2]+prob.profits[2,kp[2]]-prob.profits[2,kp[1]]]
        newW = sol.w + prob.weights[kp[2]] - prob.weights[kp[1]]
    elseif neighborhood == KP_0_1
        newX[kp[1]] = 1
        newY = [sol.y[1]+prob.profits[1,kp[1]],sol.y[2]+prob.profits[2,kp[1]]]
        newW = sol.w + prob.weights[kp[1]]
    elseif neighborhood == KP_1_2
        newX[kp[1]] = 0
        newX[kp[2]] = 1
        newX[kp[3]] = 1
        newY = [sol.y[1]+prob.profits[1,kp[2]]+prob.profits[1,kp[3]]-prob.profits[1,kp[1]],sol.y[2]+prob.profits[2,kp[2]]+prob.profits[2,kp[3]]-prob.profits[2,kp[1]]]
        newW = sol.w + prob.weights[kp[2]] + prob.weights[kp[3]] - prob.weights[kp[1]]
    end

    return Solution(newX,newY,newW)
end

function isNewSol(y::Vector{Float64}, kp::Vector{Int}, neighborhood::Neighborhood, allSolsI::Vector{Solution})
    for oldSol in allSolsI
        if oldSol.y == y
            return false
        end
    end
end

"""
    returns the index in solsInTriangles where sol should be added (triangle idex)
"""
function getTriangleIndex(solsInTriangles::Vector{Vector{Solution}}, sol::Solution)
    for i in 1:length(solsInTriangles)
        minZ1 = solsInTriangles[i][1].y[1]
        maxZ1 = solsInTriangles[i][2].y[1]
        if sol.y[1] >= minZ1 && sol.y[1] <= maxZ1
            return i,true
        end
    end
    return 0,false
    #@assert false "No triangle is fitting the solution"
end

function secureList(solsInTriangles::Vector{Vector{Solution}})
    for i in 1:length(solsInTriangles)
        for j in 3:length(solsInTriangles[i])
            for i1 in 1:length(solsInTriangles)
                for j1 in 1:length(solsInTriangles[i1])
                    if i1 != i && j != j1
                        @assert solsInTriangles[i][j].y != solsInTriangles[i1][j1].y "pas d'unicité dans la liste"
                    end
                end
            end
        end
    end
end

function whichDominates(sol1, sol2)
    sol1Dominates = true
    sol2Dominates = true

    for i in 1:2
        if sol1.y[i] < sol2.y[i]
            sol1Dominates = false
        end
        if sol1.y[i] > sol2.y[i]
            sol2Dominates = false
        end
    end

    if sol1Dominates && sol2Dominates # the Δs are the same
        @assert false "Those sols shouldn't be the same"
    elseif !sol1Dominates && !sol2Dominates # we can't compare the two kp-exchanges
        return NOTHING
    elseif sol1Dominates
        return DOMINATES
    else
        return DOMINATED
    end
end

function size2DList(list::Vector{Vector{Solution}})
    size = 0
    for i in 1:length(list)
        size += length(list[i])
    end
    return size
end

function cleanList(solsInTriangles::Vector{Vector{Solution}})
    newList = Vector{Vector{Solution}}(undef,length(solsInTriangles))
    for i in 1:length(newList)
        newList[i] = [solsInTriangles[i][1],solsInTriangles[i][2]]
    end

    fnished = false
    for i in 1:length(solsInTriangles)
        for j in 3:length(solsInTriangles[i])
            finished = false
            for i1 in 1:length(solsInTriangles)
                for j1 in 1:length(solsInTriangles[i1])
                    if i1 != i || j != j1
                        dominance_indicator = whichDominates(solsInTriangles[i][j], solsInTriangles[i1][j1])
                        if dominance_indicator == DOMINATED
                            finished = true
                            break
                        end
                        if finished
                            break
                        end
                    end
                    if finished
                        break
                    end
                end
                if finished
                    break
                end
            end

            if !finished
                push!(newList[i],solsInTriangles[i][j])
            end
        end
    end
    return newList
end

function descentNeighborhood(sol::Solution, consecutiveSet::Vector{PairOfSolution}, neighborhood::Neighborhood, prob::Problem, allSolsI::Vector{Solution}, solsInTriangles::Vector{Vector{Solution}};steepest=true,fastEvaluate=fastEvaluate)
    #println("descentNeighborhood")
    improved = false
    finished = false
    cpt = 0
    tEval = 0.
    tFun = 0.
    if neighborhood == KP_1_1
        if steepest
            bestKP,tEval,tFun = bestKP11(sol,consecutiveSet,prob,tEval,tFun,fastEvaluate=fastEvaluate)
        else
            bestKP = firstKP11(sol,consecutiveSet,prob,fastEvaluate=fastEvaluate)
        end
    elseif neighborhood == KP_0_1
        if steepest
            bestKP = bestKP01(sol,consecutiveSet,prob,fastEvaluate=fastEvaluate)
        else
            bestKP = firstKP01(sol,consecutiveSet,prob,fastEvaluate=fastEvaluate)
        end
    elseif neighborhood == KP_1_2
        if steepest
            bestKP = bestKP12(sol,consecutiveSet,prob,fastEvaluate=fastEvaluate)
        else
            bestKP = firstKP12(sol,consecutiveSet,prob,fastEvaluate=fastEvaluate)
        end
    else
        @assert false "unknown neighborhood"
    end
    push!(allSolsI, Solution(copy(sol.x),copy(sol.y),sol.w))
    while bestKP[1] != 0
        cpt += 1
        sol = updateSolByKP(sol,bestKP,neighborhood,prob)
        push!(allSolsI, Solution(copy(sol.x),copy(sol.y),sol.w))
        # test if the solution should be added to solsInTriangles and if some old ones should be deleted
        indexTriangle,triangleFound = getTriangleIndex(solsInTriangles,sol)
        if triangleFound
		    for i in 1:length(solsInTriangles[indexTriangle])
		        if solsInTriangles[indexTriangle][i].y == sol.y # the solution is not new
		            finished = true
		            break
		        end
		    end

		    if !finished
		        push!(solsInTriangles[indexTriangle],sol)
		        improved = true
		        secureList(solsInTriangles)
		    end

		    if finished
		        bestKP[1] = 0
		    elseif neighborhood == KP_1_1
		        if steepest
		            bestKP,tEval, tFun = bestKP11(sol,consecutiveSet,prob,tEval,tFun,fastEvaluate=fastEvaluate)
		        else
		            bestKP = firstKP11(sol,consecutiveSet,prob,fastEvaluate=fastEvaluate)
		        end
		    elseif neighborhood == KP_0_1
		        if steepest
		            bestKP = bestKP01(sol,consecutiveSet,prob,fastEvaluate=fastEvaluate)
		        else
		            bestKP = firstKP01(sol,consecutiveSet,prob,fastEvaluate=fastEvaluate)
		        end
		    elseif neighborhood == KP_1_2
		        if steepest
		            bestKP = bestKP12(sol,consecutiveSet,prob,fastEvaluate=fastEvaluate)
		        else
		            bestKP = firstKP12(sol,consecutiveSet,prob,fastEvaluate=fastEvaluate)
		        end
		    else
		        @assert false "Unknwon neighborhood"
		    end
	    else
	    	improved = false
	    	bestKP[1] = 0
	    end
    end

    #println("tEval = $tEval")
    #println("tFun = $tFun")

    return sol, cpt, improved
end

function descent(sol::Solution, consecutiveSet::Vector{PairOfSolution}, neighborhoods::Vector{Neighborhood}, prob::Problem, allSolsI::Vector{Solution}, solsInTriangles::Vector{Vector{Solution}}; steepest = steepest, fastEvaluate = fastEvaluate)
    #println("descent")
    cptGlobal = zeros(Int,length(neighborhoods))
    improving = false
    for i in 1:length(neighborhoods)
        sol, cpt, improved = descentNeighborhood(sol,consecutiveSet,neighborhoods[i],prob,allSolsI,solsInTriangles,steepest=steepest, fastEvaluate=fastEvaluate)
        cptGlobal[i] += cpt
        improving = improving || improved
    end

    #println("   <cptGlobal = $cptGlobal>")

    return sol, improving
end

function improveLowerBoundSet(lowerBoundSet::Vector{Solution},consecutiveSet::Vector{PairOfSolution},prob::Problem;steepest = true, fastEvaluate = true, neighborhoods = [KP_1_1])
    #println("steepest = $steepest, fastEval = $fastEvaluate, neighborhoods = $neighborhoods")
    if length(lowerBoundSet) == 0
        return lowerBoundSet, consecutiveSet
    else
        @assert length(lowerBoundSet) == length(consecutiveSet)+1 "length(lowerBoundSet) != length(consecutiveSet)+1"

        solsPrime = Vector{Solution}(undef,length(lowerBoundSet))
        allSols = Vector{Vector{Solution}}(undef,length(lowerBoundSet))
        solsInTriangles = Vector{Vector{Solution}}(undef,length(consecutiveSet))

        for i in 1:length(consecutiveSet)
            solsInTriangles[i] = [consecutiveSet[i].solL, consecutiveSet[i].solR]
        end

        for i in 1:length(lowerBoundSet)
            allSols[i] = Vector{Solution}()
        end

        improving = true
        while improving
            improving = false
            for i in 1:length(lowerBoundSet)
                solsPrime[i], improved = descent(lowerBoundSet[i],consecutiveSet, neighborhoods, prob, allSols[i], solsInTriangles, steepest = steepest, fastEvaluate = fastEvaluate)
                improving = improving || improved
            end

            # setting the structures for a new iteration
            solsClean = cleanList(solsInTriangles)
            newLowerBoundSet = Vector{Solution}()

            for i in 1:length(solsClean)
                push!(newLowerBoundSet,solsClean[i][1])
                for j in 3:length(solsClean[i])
                    push!(newLowerBoundSet,solsClean[i][j])
                end
                if i == length(solsClean)
                    push!(newLowerBoundSet,solsClean[i][2])
                end
            end

            sort!(newLowerBoundSet, by = sol -> sol.y[1])

            newConsecutiveSet = Vector{PairOfSolution}(undef,length(newLowerBoundSet)-1)
            for i in 1:length(consecutiveSet)
                consecutiveSet[i] = PairOfSolution(newLowerBoundSet[i],newLowerBoundSet[i+1])
            end
            # at this point we have constructed our new structures

        end

        #println("allSols = $allSols")

        #println("solsPrime = $solsPrime")

        #println("solsInTriangles = $solsInTriangles")
        #println("size2DList(solsInTriangles) = $(size2DList(solsInTriangles))")

        solsClean = cleanList(solsInTriangles)
        lb = nil(Solution)
        lb = cons(solsClean[length(solsClean)][2],lb)

        #=for i in 1:length(solsClean)
            print(solsClean[i])
        end=#

        lbVect = Vector{Solution}()

        for i in 1:length(solsClean)
            push!(lbVect,solsClean[i][1])
            for j in 3:length(solsClean[i])
                push!(lbVect,solsClean[i][j])
            end
            if i == length(solsClean)
                push!(lbVect,solsClean[i][2])
            end
        end

        println("lbVect = $lbVect")

        sort!(lbVect, by = sol -> sol.y[1])

        println("lbVect = $lbVect")

        lb = nil(Solution)
        lb = cons(lbVect[length(lbVect)],lb)

        for i in length(lbVect)-1:-1:1
            lb = 	cons(lbVect[i],lb)
        end

        nadirPoints = getNadirPoints(lb) #Jules

        println("nadirPoints = $nadirPoints")

        #println("solsClean = $solsClean")
        #println("length(lowerBoundSet) = $(length(lowerBoundSet))")
        #println("size2DList(solsClean) = $(size2DList(solsClean) - length(lowerBoundSet) + 2)")

        #plotRunBis("Inst200",lowerBoundSet,solsPrime, allSols, consecutiveSet)

        #plotRun3("Inst200",lowerBoundSet,solsPrime, allSols, consecutiveSet, solsClean)

        return lb, nadirPoints, size2DList(solsClean) - 2*length(lowerBoundSet) + 2, solsInTriangles
    end
end

function plotRun3(iname, sols, newSols, allSols, consecutiveSet, solsInTriangle)
    figure("Run Steepest",figsize=(6,6)) # Create a new figure
    title("Steepest-Descent-2OKP | " * iname)
    xlabel("z1")
    ylabel("z2")

    #println("consecutiveSet = $consecutiveSet")
    for i in 1:length(solsInTriangle)
        for j in 3:length(solsInTriangle[i])
            sol = solsInTriangle[i][j]

            plot(sol.y[1], sol.y[2], ms = "4", marker = "o", color = "black")
        end
    end
end

function plotRunBis(iname, sols, newSols, allSols, consecutiveSet)
    figure("Run Steepest",figsize=(6,6)) # Create a new figure
    #title("Steepest-Descent-2OKP | " * iname)
    xlabel("z1")
    ylabel("z2")

    x = Vector{Float64}(undef,length(sols))
    y = Vector{Float64}(undef,length(sols))
    newX = Vector{Float64}(undef,length(newSols))
    newY = Vector{Float64}(undef,length(newSols))
    for i in 1:length(sols)
        x[i] = sols[i].y[1]
        y[i] = sols[i].y[2]
    end

    for i in 1:length(newSols)
        newX[i] = newSols[i].y[1]
        newY[i] = newSols[i].y[2]
    end

    X = Vector{Vector{Float64}}(undef,length(allSols))
    Y = Vector{Vector{Float64}}(undef,length(allSols))
    for i in 1:length(allSols)
        X[i] = Vector{Float64}()
        Y[i] = Vector{Float64}()
        for j in 1:length(allSols[i])
            push!(X[i],allSols[i][j].y[1])
            push!(Y[i],allSols[i][j].y[2])
        end
    end

    XT = Vector{Float64}(undef,3*length(consecutiveSet))
    YT = Vector{Float64}(undef,3*length(consecutiveSet))
    #println("consecutiveSet = $consecutiveSet")
    for i in 1:length(consecutiveSet)
        pair = consecutiveSet[i]
        solL = pair.solL
        solR = pair.solR

        plot([solL.y[1], solL.y[1], solR.y[1], solL.y[1]], [solL.y[2], solR.y[2], solR.y[2], solL.y[2]], color = "orange", ls = "-", linewidth=1)

        XT[3*(i-1)+1] = solL.y[1]
        XT[3*(i-1)+2] = solL.y[1]
        XT[3*(i-1)+3] = solR.y[1]

        YT[3*(i-1)+1] = solL.y[2]
        YT[3*(i-1)+2] = solR.y[2]
        YT[3*(i-1)+3] = solR.y[2]
    end

    for i in 1:length(allSols)
        plot(X[i],Y[i],ls="-", marker = "o", ms = 3)
    end

    #plot(XT,YT,ls="-",marker="o",ms=3,color="blue",label="Triangles")
    plot(x,y,ls="",marker="o",ms=4,color="orange",label="LB")
    #plot(newX,newY,ls="",marker="x",ms=10,color="red",label="nouvelles")
    #legend(loc=4, fontsize ="small")
end

function plotRun(iname, sols, newSols)
    figure("Run",figsize=(6,6)) # Create a new figure
    #title("Steepest-Descent-2OKP | " * iname)
    xlabel("z1")
    ylabel("z2")

    x = Vector{Float64}(undef,length(sols))
    y = Vector{Float64}(undef,length(sols))
    newX = Vector{Float64}(undef,length(newSols))
    newY = Vector{Float64}(undef,length(newSols))
    for i in 1:length(sols)
        x[i] = sols[i].y[1]
        y[i] = sols[i].y[2]
    end

    for i in 1:length(newSols)
        newX[i] = newSols[i].y[1]
        newY[i] = newSols[i].y[2]
    end

    #plot(x,y,ls="",marker="+",ms=5,color="green",label="LB")
    #plot(newX,newY,ls="",marker="x",ms=5,color="red",label="nouvelles")
    #legend(loc=4, fontsize ="small")
end

max_stars_repo_name: LucasBaussay/BranchAndBoundBiObjKp

max_issues_repo_name: LucasBaussay/BranchAndBoundBiObjKp

max_forks_repo_name: LucasBaussay/BranchAndBoundBiObjKp


text:
"""Running scenarios one after the other."""

using DataFrames, CSV, FileIO
using BenchmarkTools
using Agtor, Dates

import Base.Filesystem: rm

using Statistics
using PyCall

# Start month/day of growing season
# Represents when water allocations start getting announced.
const gs_start = (5, 1)

function setup_basin(data_dir, climate_data)
    basin_spec_dir::String = "$(data_dir)basins"
    basin_spec::Dict = load_spec(basin_spec_dir)[:TestBasin]

    basin_name = basin_spec[:name]
    zone_specs = basin_spec[:zone_spec]

    Manager_A = BaseManager("optimizing")
    Manager_B = RigidManager("rigid", 0.05)
    manage_zones = ((Manager_A, ("Zone_1", )), (Manager_B, ("Zone_2", )))

    b = Basin(name=basin_name, zone_spec=zone_specs, climate_data=climate_data, managers=manage_zones)

    return b
end

function calculate_metrics(collated_results::Union{Tuple, Array})::Array{Float64}

    z_res = DataFrame(fill(Float64, 4), Symbol.((i for i in range(1,4))), length(collated_results))

    for (idx, zone_res) in enumerate(collated_results)
        z_res[idx, :] = calculate_metrics(zone_res)
    end

    return [mean(z_res[:, i]) for i in names(z_res)]
end

function calculate_metrics(collated_logs::DataFrame)
    income = collated_logs[!, Symbol("income_sum")]

    cv = std(income) / mean(income)
    mean_ML_per_yield = mean(collated_logs[!, Symbol("ML per Irrigated Yield")])
    mean_income_per_ha = mean(collated_logs[!, Symbol("Dollar per Ha")])
    mean_water_use = mean(collated_logs[!, :irrigated_volume_sum])

    return [cv, mean_ML_per_yield, mean_income_per_ha, mean_water_use]
end

function run_model(scenario_id, basin)

    # Example annual water allocations
    allocs = (surface_water=150.0, groundwater=40.0)

    for dt_i in basin.climate.time_steps
        # Resetting allocations for example run
        if monthday(dt_i) == gs_start
            for z in basin.zones
                update_available_water!(z, allocs)
            end
        end

        run_timestep!(basin)
    end

    # Save results as they complete
    all_res = []
    for (z_id, z) in enumerate(basin.zones)
        results = collect_results(z)

        # run_id = "$(scenario_id)/$(z.name)"
        # save_results!("basin_async_save.jld2", run_id, results; mode="a+")

        push!(all_res, results[1])
    end

    return calculate_metrics(all_res)
end

function scenario_run(samples, basin)::Array{Array{Float64}}
    basin_wide_res = []
    for (row_id, r) in enumerate(eachrow(samples))
        tmp_b = deepcopy(basin)

        update_model!(tmp_b, r)
        push!(basin_wide_res, run_model(row_id, tmp_b))

        # An alternative to using a deepcopy of the basin is
        # to reset its values (more flexible, but increases runtime slightly)
        # reset!(basin)
    end

    return basin_wide_res
end

data_dir = "test/data/"
example_basin = setup_basin(data_dir, "$(data_dir)climate/basin_climate_data.csv")

all_params = collect_agparams!(example_basin, [])
param_names = all_params[:name]
param_bounds = Array(all_params[:, [:min_val, :max_val]])

py"""
from SALib import ProblemSpec
import numpy as np
import pandas as pd

sp = ProblemSpec({
  'names': $(param_names),
  'bounds': $(Array(param_bounds)),
  'output': ['CV', 'Mean ML per Yield [ML/t]', 'Mean Income [$]', 'Mean Water Use [ML]']
})

sp.sample_latin(4, seed=101)
samples = pd.DataFrame(sp.samples, columns=sp['names'])
samples.to_csv("basin_samples.csv", index=False)
"""

samples = DataFrame!(CSV.File("basin_samples.csv"))
@time res = scenario_run(samples, example_basin)

py"""
sp.results = np.array($(res))
sp.analyze_rbd_fast()
"""

@info "Analysis:" py"sp.to_df()"

rm("basin_samples.csv")

ext: jl

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
@model school8(J, y, sigma) = begin
  mu ~ Flat()
  tau ~ FlatPos(0)
  eta = tzeros(J)
  eta ~ [Normal(0, 1)]
  y ~ MvNormal(mu .+ tau .* eta, sigma)
end


text:
module Bayesian_Survival_Analysis

include("generate_dummy_data.jl")

export generateold, generatecensored, generategakui, generategender, generateseniority, generateyearofjoin

end # module


text:
# tests for problems and ensemble problem types

@testset "StateTransfer in place w/ GRAPE" begin

    prob = Problem(
        B = [Sx, Sy],
        A = Sz,
        Xi = ρinit,
        Xt = ρfin,
        T = 1.0,
        n_controls = 2,
        guess = rand(2, 10),
        sys_type = StateTransfer(),
    )

    sol = solve(prob, GRAPE(n_slices = 10, isinplace = true))

    @test sol.result.minimum - C1(ρfin, ρfin) < tol
end

@testset "StateTransfer GRAPE" begin

    prob = Problem(
        B = [SSx, SSy],
        A = SSz,
        Xi = ρinitS,
        Xt = ρfinS,
        T = 1.0,
        n_controls = 2,
        guess = rand(2, 10),
        sys_type = StateTransfer(),
    )

    sol = solve(prob, GRAPE(n_slices = 10, isinplace = false))
    @test sol.result.minimum - C1(ρfinS, ρfinS) < tol

end

@testset "StateTransfer Ensemble inplace" begin

    prob = Problem(
        B = [Sx, Sy],
        A = Sz,
        Xi = ρinit,
        Xt = ρfin,
        T = 5.0,
        n_controls = 2,
        guess = rand(2, 25),
        sys_type = StateTransfer(),
    )

    ens = EnsembleProblem(
        prob = prob,
        n_ens = 5,
        A_g = A_gens,
        B_g = B_gens,
        XiG = X_init_gens,
        XtG = X_target_gens,
        wts = ones(5) / 5,
    )

    sol = solve(ens, GRAPE(n_slices = 25, isinplace = true))
    @test sol.result.minimum - C1(ρfin, ρfin) < tol * 10

end

@testset "StateTransfer Ensemble out of place" begin

    prob = Problem(
        B = [SSx, SSy],
        A = SSz,
        Xi = ρinitS,
        Xt = ρfinS,
        T = 5.0,
        n_controls = 2,
        guess = rand(2, 25),
        sys_type = StateTransfer(),
    )

    ens = EnsembleProblem(
        prob = prob,
        n_ens = 5,
        A_g = A_gens_static,
        B_g = B_gens_static,
        XiG = X_init_gens_static,
        XtG = X_target_gens_static,
        wts = ones(5) / 5,
    )

    sol = solve(ens, GRAPE(n_slices = 25, isinplace = false))
    @test sol.result.minimum - C1(ρfinS, ρfinS) < tol * 10

end

@testset "StateTransfer w/ ADGRAPE" begin

    prob = Problem(
        B = [Sx, Sy],
        A = Sz,
        Xi = ρinit,
        Xt = ρfin,
        T = 1.0,
        n_controls = 2,
        guess = rand(2, 10),
        sys_type = StateTransfer(),
    )

    sol = solve(prob, ADGRAPE(n_slices = 10))

    @test sol.result.minimum - C1(ρfin, ρfin) < tol
end

@testset "StateTransfer Ensemble ADGRAPE" begin

    prob = Problem(
        B = [Sx, Sy],
        A = Sz,
        Xi = ρinit,
        Xt = ρfin,
        T = 5.0,
        n_controls = 2,
        guess = rand(2, 25),
        sys_type = StateTransfer(),
    )

    ens = EnsembleProblem(
        prob = prob,
        n_ens = 5,
        A_g = A_gens,
        B_g = B_gens,
        XiG = X_init_gens,
        XtG = X_target_gens,
        wts = ones(5) / 5,
    )

    sol = solve(ens, ADGRAPE(n_slices = 25))
    @test sol.result.minimum - C1(ρfin, ρfin) < tol

end


text:
#   Generating #define constants
const XML_TRUE = 1
const XML_FALSE = 0


text:
module TestIndex

using Test, DataFrames
using DataFrames: Index, SubIndex, fuzzymatch

i = Index()
push!(i, :A)
push!(i, :B)

inds = Any[1,
           big(1),
           :A,
           [true, false],
           [1],
           [big(1)],
           big(1):big(1),
           [:A],
           Union{Int, Missing}[1],
           Union{BigInt, Missing}[big(1)],
           Union{Symbol, Missing}[:A],
           Any[1],
           Any[:A]]

for ind in inds
    if ind == :A || ndims(ind) == 0
        @test i[ind] == 1
    else
        @test (i[ind] == [1])
    end
end

@test_throws MethodError i[1.0]
@test_throws ArgumentError i[true]
@test_throws ArgumentError i[false]
@test_throws ArgumentError i[Union{Bool, Missing}[true, false]]
@test_throws ArgumentError i[Any[1, missing]]
@test_throws ArgumentError i[[1, missing]]
@test_throws ArgumentError i[[true, missing]]
@test_throws ArgumentError i[Any[true, missing]]
@test_throws MethodError i[[:A, missing]]
@test_throws MethodError i[Any[:A, missing]]
@test_throws ArgumentError i[1.0:1.0]
@test_throws ArgumentError i[[1.0]]
@test_throws ArgumentError i[Any[1.0]]

@test i[1:1] == 1:1

@test_throws BoundsError i[[true]]
@test_throws BoundsError i[true:true]
@test_throws BoundsError i[[true, false, true]]

@test_throws ArgumentError i[["a"]]
@test_throws ArgumentError i[Any["a"]]

@test i[[]] == Int[]
@test i[Int[]] == Int[]
@test i[Symbol[]] == Int[]

@test names(i) == [:A,:B]
@test names!(i, [:a,:a], makeunique=true) == Index([:a,:a_1])
@test_throws ArgumentError names!(i, [:a,:a])
@test names!(i, [:a,:b]) == Index([:a,:b])
@test rename(i, Dict(:a=>:A, :b=>:B)) == Index([:A,:B])
@test rename(i, :a => :A) == Index([:A,:b])
@test rename(i, :a => :a) == Index([:a,:b])
@test rename(i, [:a => :A]) == Index([:A,:b])
@test rename(i, [:a => :a]) == Index([:a,:b])
@test rename(x->Symbol(uppercase(string(x))), i) == Index([:A,:B])
@test rename(x->Symbol(lowercase(string(x))), i) == Index([:a,:b])

@test delete!(i, :a) == Index([:b])
push!(i, :C)
@test delete!(i, 1) == Index([:C])

i = Index([:A, :B, :C, :D, :E])
i2 = copy(i)
names!(i2, reverse(names(i2)))
names!(i2, reverse(names(i2)))
@test names(i2) == names(i)
for name in names(i)
  i2[name] # Issue #715
end

i = Index([:A, :B, :C, :D, :E])
si1 = SubIndex(i, :)
si2 = SubIndex(i, 3:5)
si3 = SubIndex(i, [3,4,5])
si4 = SubIndex(i, [false, false, true, true, true])
si5 = SubIndex(i, [:C, :D, :E])

@test copy(si1) == i
@test copy(si2) == Index([:C, :D, :E])
@test copy(si3) == Index([:C, :D, :E])
@test copy(si4) == Index([:C, :D, :E])
@test copy(si5) == Index([:C, :D, :E])

@test_throws ArgumentError SubIndex(i, 1)
@test_throws ArgumentError SubIndex(i, :A)
@test_throws ArgumentError SubIndex(i, true)
@test si1 isa Index
@test si2.cols == 3:5
@test si2.remap == -1:3
@test si3.cols == 3:5
@test si3.remap == Int[]
@test !haskey(si3, :A)
@test si3.remap == [0, 0, 1, 2, 3]
@test si4.cols == 3:5
@test si4.remap == Int[]
@test !haskey(si4, :A)
@test si4.remap == [0, 0, 1, 2, 3]
@test si5.cols == 3:5
@test si5.remap == Int[]
@test !haskey(si5, :A)
@test si5.remap == [0, 0, 1, 2, 3]

@test length(si1) == 5
@test length(si2) == 3
@test length(si3) == 3
@test length(si4) == 3
@test length(si5) == 3

@test names(si1) == keys(si1) == [:A, :B, :C, :D, :E]
@test names(si2) == keys(si2) == [:C, :D, :E]
@test names(si3) == keys(si3) == [:C, :D, :E]
@test names(si4) == keys(si4) == [:C, :D, :E]
@test names(si5) == keys(si5) == [:C, :D, :E]

@test_throws ArgumentError haskey(si3, true)
@test haskey(si3, 1)
@test !haskey(si3, 0)
@test !haskey(si3, 4)
@test haskey(si3, :D)
@test !haskey(si3, :A)
@test si3[:C] == 1
@test si3[names(i)] == [0, 0, 1, 2, 3]

@testset "selector mutation" begin
    df = DataFrame(a=1:5, b=11:15, c=21:25)
    selector1 = [3,2]
    dfv1 = view(df, selector1)
    dfr1 = view(df, 2, selector1)
    selector2 = [1]
    dfv2 = view(dfv1, selector2)
    dfr2 = view(dfr1, selector2)
    @test names(dfv1) == [:c, :b]
    @test names(dfv2) == [:c]
    @test names(dfr1) == [:c, :b]
    @test names(dfr2) == [:c]
    selector1[1] = 1
    @test names(dfv1) == [:a, :b]
    @test names(dfv2) == [:c]
    @test names(dfr1) == [:a, :b]
    @test names(dfr2) == [:c]
    selector3 = [:c, :b]
    dfv3 = view(df, selector3)
    dfr3 = view(df, 2, selector3)
    @test names(dfv3) == [:c, :b]
    @test names(dfr3) == [:c, :b]
    selector3[1] = :a
    @test names(dfv3) == [:c, :b]
    @test names(dfr3) == [:c, :b]
end

@testset "fuzzy matching and ArgumentError" begin
    i = Index()
    push!(i, :x1)
    push!(i, :x12)
    push!(i, :x131)
    push!(i, :y13)
    push!(i, :yy13)
    @test_throws ArgumentError i[:x13]
    @test_throws ArgumentError i[:xx13]
    @test_throws ArgumentError i[:yy14]
    @test_throws ArgumentError i[:abcd]
    @test fuzzymatch(i.lookup, :x13) == [:x1, :x12, :x131, :y13, :yy13]
    @test fuzzymatch(i.lookup, :xx1314) == [:x131]
    @test fuzzymatch(i.lookup, :yy14) == [:yy13, :y13]
    @test isempty(fuzzymatch(i.lookup, :abcd))
end

end # module


text:
"""

units (x,y,z) are in millimeters; z assumed vertical (column length)
density in thousands per cubic millimeter
λ -- connection probability decay rate
v -- conduction velocity (m/s (or mm/ms))
"""

function analyzeblock(delmat::Array{Int64,2}, fs)
	n = size(delmat)[1]
	conmat = map(x -> x != 0 ? 1 : 0, delmat)
	numdels = sum(conmat)
	println("Number of nodes:\t$(n)")
	println("Probability of contact:\t$(numdels/(n*n))")
	println("Average delay:\t$(1000.0*sum(delmat)/(fs*numdels)) ms")
end

function analyzeblock(delmat::Array{Int64,2}, fs, ps,
					  bounds1::Tuple{Float64,Float64,Float64},
					  bounds2::Tuple{Float64,Float64,Float64})
	n = size(delmat)[1]
	idcs = findall( x -> bounds1[1] <= x[1] <= bounds2[1] &&
				   		 bounds1[2] <= x[2] <= bounds2[2] &&
				   		 bounds1[3] <= x[3] <= bounds2[3],
						 ps)
	conmat = map(x -> x != 0 ? 1 : 0, delmat)
	numdels = sum(conmat[idcs,:])
	p_contact = numdels/(length(idcs)*length(idcs))
	avgdelay = 1000.0*sum(delmat[idcs,:])/(fs*numdels)
	println("Number of nodes:\t$(length(idcs))")
	println("Probability of contact:\t$p_contact")
	println("Average delay:\t$avgdelay ms")
end

function delayhistvals(delmat::Array{Int64,2}, fs)
	n = size(delmat)[1]
	numdels = sum( map( x -> x != 0 ? 1 : 0, delmat ) )
	vals = zeros(Float64, numdels)
	idx = 1
	for i ∈ 1:size(delmat)[1]
		for j ∈ 1:size(delmat)[2]
			if delmat[i,j] != 0
				vals[idx] = 1000.0 * delmat[i,j] / fs
				idx += 1
			end
		end
	end
	vals
end

function lenstobounds(lens::Array{Int64,1})
	bounds = Array{Int64,1}(undef,length(lens)+1)
	bounds[1] = 1
	for i ∈ 1:length(lens)
		bounds[i+1] = bounds[i] + lens[i]
	end
	bounds
end

dist(p1, p2) = sqrt(sum((p1 .- p2) .^ 2))
delaysamples(d,v,fs) = Int(round((d/v)*fs))

function sampletargets(center, positions, λ; factor=1.0) # λ = scale length
    n = length(positions)
    targets = zeros(Int64,0)
    for i ∈ 1:n
        d = dist(center, positions[i])
	    p_contact = factor*λ*exp( -d/λ )
        if rand() < p_contact
            push!(targets, i)
        end
    end
    targets
end

function sampletargets(center, positions, λ,
                       neurontypes, type::SimpleNeuronType;
                       factor=1.0)
    n = length(positions)
    targets = zeros(Int64,0)
    for i ∈ 1:n
        d = dist(center, positions[i])
	    p_contact = factor*λ*exp( -d/λ )
        if (rand() < p_contact) && (neurontypes[i] == type)
            push!(targets, i)
        end
    end
    targets
end

function genpatch(dims::Tuple{T,T,T},
				  types::Array{SimpleNeuronType,1},
				  ρs::Array{Float64,1},
				  λs::Array{Float64,1},
				  vs::Array{Float64,1},
				  fs;
				  probfactor = [1.0 for _ ∈ 1:length(types)],
				  numslices = 10,
				  maxlen=1.0,
				  verbose=false) where T <: Number
	@assert length(types) == length(ρs) == length(λs) == length(vs)
	slicedepth = dims[3]/numslices;
	positions = Array{Array{Tuple{T,T,T},1},1}(undef,0)
	neurons = Array{Array{SimpleNeuronType,1},1}(undef,0)
	lens = Array{Int64,1}(undef,length(types))
	for (i,ρ) ∈ enumerate(ρs)
		pointsinslice = dims[1]*dims[2]*slicedepth*ρs[i]
		ps = vcat([[(rand()*dims[1], rand()*dims[2], slicedepth*(rand()+k))
					for _ ∈ 1:pointsinslice] for k ∈ 0:numslices-1]...)
		push!(positions, ps)
		ns = [types[i] for _ ∈ 1:length(ps)]
		push!(neurons,ns)
		lens[i] = length(ns)
	end
	positions = vcat(positions...)
	neurontypes = vcat(neurons...)
	n = length(positions)
	delgraph = zeros(Int64,n,n)
	idx = 1
	bounds = lenstobounds(lens)
	for i ∈ 1:n
		if i == bounds[idx+1]
			idx += 1
		end
		for j ∈ 1:n
			d = dist(positions[i], positions[j])
			p_contact = probfactor[idx]*exp( -d/λs[idx] )
			if d < maxlen && rand() < p_contact
				if neurontypes[i] == SimpleNeuronType(:fs) ||
					neurontypes[i] == SimpleNeuronType(:inh)
					if neurontypes[i] == neurontypes[j]
						delgraph[i,j] = 0 #no inh->inh connections
					else
						#delgraph[i,j] = 1
						delgraph[i,j] = max(delaysamples(d,vs[idx],fs), 1)
					end
				else
					delgraph[i,j] = max(delaysamples(d,vs[idx],fs), 1)
					delgraph[i,j] = delgraph[i,j] == 1 ? 1 + Int(round(rand()*4.0)) : delgraph[i,j]
				end
			end
		end
	end

	verbose && (analyzeblock(delgraph,fs))
	verbose && (println("Average distance (overall): ", avgdistance(positions)))
	verbose && (println("Average distance (connected): ",
						avgdistance(positions, delgraph)))

	return positions, neurontypes, delgraph
end

function gensyn(delgraph, neurontypes, weights::Dict{SimpleNeuronType,Float64})
	@assert length(neurontypes) == size(delgraph)[1]
	syngraph = zeros(Float64, size(delgraph))
	for i ∈ 1:size(delgraph)[1]
		w = weights[neurontypes[i]]
		for j ∈ 1:size(delgraph)[2]
			syngraph[i,j] = delgraph[i,j] != 0 ? w : 0.0
		end
	end
	syngraph
end

function avgdistance(ps)
	n = length(ps)
	cumulativedist = 0.0
	for i ∈ 1:n
		for j ∈ 1:n
			cumulativedist += sqrt( sum( (ps[i] .- ps[j]) .^ 2 ) )
		end
	end
	cumulativedist / (n * n)
end

function avgdistance(ps, graph)
	n = length(ps)
	cumulativedist = 0.0
	count = 0
	for i ∈ 1:n
		for j ∈ 1:n
			if graph[i,j] != 0
				cumulativedist += sqrt( sum( (ps[i] .- ps[j]) .^ 2 ) )
				count += 1
			end
		end
	end
	cumulativedist / count
end

function distances(ps, graph)
	n = length(ps)
	distances = Array{Float64,1}(undef,n*n)
	count = 1
	for i ∈ 1:n
		for j ∈ 1:n
			if graph[i,j] != 0
				distances[count] = sqrt( sum( (ps[i] .- ps[j]) .^ 2 ) )
				count += 1
			end
		end
	end
	resize!(distances, count-1)
end

max_stars_repo_name: analogouscircuit/DelayNetTools

max_issues_repo_name: analogouscircuit/DelayNetTools

max_forks_repo_name: analogouscircuit/DelayNetTools