OrdonnancementCorrectif/src/OrdoCorr/Solver/TrackPlansILP/SimulatedAnnealing.cpp

#include "SimulatedAnnealing.hpp"
#include "../PseCarlierRivreau/omp.h"
#include "DynamicProgramming.hpp"
#include "Solution.hpp"
#include "TrackPlan.hpp"
#include <algorithm>
#include <array>
#include <cmath>
#include <ctime>
#include <iterator>
#include <optional>
#include <random>
#include <string>
#include <utility>
#include <vector>
#include "Random.hpp"
#include "sourceSolTrPlan.hpp"

namespace solverlib {
    using namespace random;

    bool StatSimulatedAnnealing::activate = true;
    bool SimulatedAnnealing::withDynProg = false;


    std::unordered_map<EMovingOperators, std::string> StatSimulatedAnnealing::names = {
        {EMovingOperators::CHANGE_MODE_WITHOUT_CARLIER, "CHANGE_MODE"},
        {EMovingOperators::INSERT_WITHOUT_CARLIER, "INSERT"},
        {EMovingOperators::REMOVE_WITHOUT_CARLIER, "REMOVE"},
        {EMovingOperators::SWAP_WITHOUT_CARLIER, "SWAP"},
        {EMovingOperators::SWAP_WITHIN_INTERVAL, "SWAP_WITHIN_SEQUENCE"},
        {EMovingOperators::DYN_PROG, "DYN_PROG"},
        {EMovingOperators::MOVE, "MOVE"}
    };


    SimulatedAnnealing::SimulatedAnnealing(std::unordered_map<unsigned short, Decision>& decs, std::shared_ptr<modellib::STFMockInstance> mock, ESourceTrackPlan source)
    {
        randomEngine = solverlib::random::makeEngine();
        addSolutionToPool(decs, mock, source);
    }

    void SimulatedAnnealing::addSolutionToPool(std::unordered_map<unsigned short, Decision>& decs, std::shared_ptr<modellib::STFMockInstance> mock, ESourceTrackPlan source, bool isPutFirst)
    {
        auto costs = evaluate(decs);
        if(!isPutFirst)
            solutions.push_back({source, mock, decs, costs.first, costs.second});
        else
        {
            solutions.insert(solutions.begin(), {source, mock, decs, costs.first, costs.second});
        }
    }


    std::pair<unsigned int, unsigned int> SimulatedAnnealing::evaluate(const std::unordered_map<unsigned short, Decision>& decs)
    {
        unsigned int cost = 0;
        unsigned int diagCost = 0;
        for(auto& dec : decs)
        {
            if(dec.second.excluded)
            {
                cost += MAXIMUM_TIME_OFFSET * modellib::STFMockInstance::jobs[dec.first]->getPoidsRetard();
            }
            else
                cost += modellib::STFMockInstance::jobs[dec.first]->getPoidsRetard() * dec.second.lastCreneau.first;

            if(dec.second.rejected)
            {
                diagCost += modellib::STFMockInstance::jobs[dec.first]->getPoidsRejet();
            }
        }
        return {cost, diagCost};
    }

    std::pair<unsigned int, unsigned int> SimulatedAnnealing::getCostOfSequence(std::vector<std::pair<unsigned short, decision>>& jobsSeq)
    {
        unsigned int cost = 0;
        unsigned int costDiag = 0;
        for(auto& el : jobsSeq)
        {
            cost += STFMockInstance::jobs[el.first]->getPoidsRetard()*el.second.lastCreneau.first;
            costDiag += STFMockInstance::jobs[el.first]->getPoidsRejet()*el.second.rejected;
        }
        return {cost, costDiag};
    }


    // tire un opérateur uniformément
    EMovingOperators SimulatedAnnealing::pick_operator(double temp, double tmax) {
        int N = withDynProg ? static_cast<int>(EMovingOperators::DYN_PROG)+1 : static_cast<int>(EMovingOperators::MOVE)+1;
        std::vector<float> weights(N, 1);

        std::vector<float> base_weights = {
            1.0f,   
            1.0f,   
            1.0f,
            1.0f,   
            1.0f,   
            1.0f    
            // DYN_PROG ajouté si besoin
        };
        if (withDynProg) base_weights.push_back(0.5f); // DYN_PROG

        const float remove_base = 1.0f;  // poids max du remove (en début de recuit)
        const float t = temp / tmax;  // 1.0 → 0.0

        // Poids remove
        const float w_remove = remove_base * t;

        // Le budget récupéré est redistribué proportionnellement aux autres
        const float base_sum = std::accumulate(base_weights.begin(), base_weights.end(), 0.0f);
        const float bonus = remove_base * (1.0f - t);

        // Index du REMOVE dans ton enum — à adapter
        constexpr int REMOVE_IDX = 2;

        for (int i = 0; i < N; ++i) {
            if (i == REMOVE_IDX) {
                weights[i] = w_remove;
            } else {
                weights[i] = base_weights[i] + bonus * (base_weights[i] / base_sum);
            }
        }
        return static_cast<EMovingOperators>(
            //std::uniform_int_distribution<int>(0, N)(randomEngine)
            std::discrete_distribution<int>(weights.begin(), weights.end())(randomEngine)
        );
    }

    // Applique un opérateur et retourne un voisin (nullopt si infaisable)
    std::optional<SASolution> SimulatedAnnealing::apply_operator(const SASolution& current, EMovingOperators op)
    {
        switch (op)
        {
            case EMovingOperators::SWAP_WITHOUT_CARLIER:                    return move_swap_WC(current);
            case EMovingOperators::INSERT_WITHOUT_CARLIER:                  return move_insert_WC(current);
            case EMovingOperators::REMOVE_WITHOUT_CARLIER:                  return move_remove_WC(current);
            case EMovingOperators::CHANGE_MODE_WITHOUT_CARLIER:             return move_change_mode_WC(current);
            case EMovingOperators::SWAP_WITHIN_INTERVAL:                    return move_swap_within_interval(current);
            case EMovingOperators::MOVE:                                    return move_move_WC(current);
            case EMovingOperators::DYN_PROG:                                return move_dynprog(current);                       
              break;
            }
        return std::nullopt;
    }

    // Boucle principale
    SASolution SimulatedAnnealing::solve(double T_max, double T_min, double cooling_rate, int iterations_per_temp)
    {
        SASolution current  = solutions[0];
        SASolution best     = solutions[0];
        double   cost_cur = current.cost;
        double   cost_best= cost_cur;
        double   T        = T_max;
        std::uniform_real_distribution<double> uniform(0.0, 1.0);

        while (T > T_min + 10e-6)
        {
            for (int i = 0; i < iterations_per_temp; ++i)
            {
                EMovingOperators op = pick_operator(T, T_max);
                stats.addUsed(op);
                auto neighbor = apply_operator(current, op);

                if (!neighbor.has_value()) continue;

                auto costs_neighbor = std::make_pair(neighbor->cost, neighbor->diagCost);

                double delta = costs_neighbor.first - cost_cur;

                stats.addFeas(op);
                if(delta <= 0)
                    stats.addImproved(op);

                if(delta > 0)
                {
                    stats.addFailInfo(op, getP(delta, T, op), delta, T);   
                }

                if (delta < 0 || uniform(randomEngine) < getP(delta, T, op))
                {
                    current  = std::move(*neighbor);
                    cost_cur = costs_neighbor.first;

                    if (cost_cur < cost_best) {
                        best      = current;
                        cost_best = cost_cur;
                    }

                    solutions.push_back(current);
                }
            }

            T *= cooling_rate;
            //loggerlib::Logger::systemNotify(loggerlib::LOGGER_PROGRESS, "Temperature : " + std::to_string(T));
            //loggerlib::Logger::systemNotify(loggerlib::LOGGER_PROGRESS, "Solution pool : " + std::to_string(solutions.size()));
        }
        //loggerlib::Logger::systemNotify(loggerlib::LOGGER_PROGRESS, "Best solution : " + std::to_string(best.cost));

        return best;
    }

    double SimulatedAnnealing::getP(double delta, double temperature, EMovingOperators op)
    {
        switch (op) {
            case EMovingOperators::SWAP_WITHOUT_CARLIER:
            case EMovingOperators::INSERT_WITHOUT_CARLIER:
            case EMovingOperators::REMOVE_WITHOUT_CARLIER: return std::exp(-delta/(temperature*100));
            case EMovingOperators::CHANGE_MODE_WITHOUT_CARLIER:
            case EMovingOperators::SWAP_WITHIN_INTERVAL:
            case EMovingOperators::MOVE: return std::exp(-delta/(temperature*10));
            case EMovingOperators::DYN_PROG:
          break;
        }
        return std::exp(-delta/temperature);
    }


    //OPERATEURS

    std::optional<SASolution> SimulatedAnnealing::move_dynprog(const SASolution& sol)
    {
        DynamicProgramming prog(sol);
        prog.mode = true;
        prog.saveSols = false;
        auto result = prog.solve();
        return std::optional<SASolution>(result);
    }

    std::optional<SASolution> SimulatedAnnealing::move_swap_within_interval(const SASolution& sol)
    {
        auto mock = sol.mock;
        if (!mock) return std::nullopt;

        std::vector<unsigned short> active_ops;
        for (auto& [op_id, dec] : sol.decisions)
            if (!dec.excluded) active_ops.push_back(op_id);
        if (active_ops.size() < 2) return std::nullopt;

        std::uniform_int_distribution<int> dist(0, (int)active_ops.size() - 1);
        unsigned short op_a = active_ops[dist(randomEngine)];
        const Decision& dec_a = sol.decisions.at(op_a);

        std::vector<unsigned short> seq_swap(1, op_a);

        for (auto& op_id : active_ops) {
            if (op_id == op_a) continue;
            const Decision& dec_b = sol.decisions.at(op_id);
            if (dec_b.empV != dec_a.empV) continue;

            seq_swap.push_back({
                op_id
            });
        }
        if (seq_swap.size() == 1) return std::nullopt;

        std::sort(seq_swap.begin(), seq_swap.end(), [&](auto job1, auto job2){
            const Decision& dec_a = sol.decisions.at(job1);
            const Decision& dec_b = sol.decisions.at(job2);

            return dec_a.lastCreneau.first < dec_b.lastCreneau.first;
        });

        auto posA = std::find(seq_swap.begin(), seq_swap.end(), op_a);
        if(posA == seq_swap.end()) return std::nullopt;
        
        auto lbdBuildSeq = [&](std::vector<unsigned short>& seq) -> std::optional<SASolution> {
            SASolution neighbor = sol;
            auto creneauTrack = STFMockInstance::machines[dec_a.empV]->getDispo();
            unsigned short lastEnd = 0;
            unsigned int oldCost = 0;
            unsigned int newCost = 0;
            for(auto& job : seq)
            {
                Decision& dec_neih = neighbor.decisions.at(job);
                auto creneauJob = mock->trajectoryStops[dec_neih.empR].getDispoStop();
                auto match = CreneauHoraire::checkSlotsCompatibility(creneauJob, creneauTrack);
                oldCost += STFMockInstance::jobs[job]->getPoidsRetard()*sol.decisions.at(job).lastCreneau.first;
                if(match.first)
                {
                    if(std::max(match.second.first, lastEnd) + dec_neih.rejected*STFMockInstance::jobs[job]->getDureeDiag() + !dec_neih.rejected*STFMockInstance::jobs[job]->getDuree() > match.second.first + match.second.second)
                        return std::nullopt;

                    dec_neih.lastCreneau = {
                        std::max(match.second.first, lastEnd),
                        std::max(match.second.first, lastEnd)+ dec_neih.rejected*STFMockInstance::jobs[job]->getDureeDiag() + !dec_neih.rejected*STFMockInstance::jobs[job]->getDuree()
                    };
                    newCost += STFMockInstance::jobs[job]->getPoidsRetard()*dec_neih.lastCreneau.first;
                    lastEnd = std::max(match.second.first, lastEnd)+ dec_neih.rejected*STFMockInstance::jobs[job]->getDureeDiag() + !dec_neih.rejected*STFMockInstance::jobs[job]->getDuree();
                }
                else {
                    return std::nullopt;
                }
            }
            neighbor.cost = (neighbor.cost - oldCost) + newCost;
            neighbor.source = ESourceTrackPlan::SimAn;
            return neighbor;
        };

        //TEST SWAP
        std::uniform_int_distribution<int> posR(0, (int)seq_swap.size()-1);
        auto job = posR(randomEngine);
        auto IposA = std::distance(seq_swap.begin(), posA);
        while(job == IposA)
            job = posR(randomEngine);

        auto seqCop = seq_swap;
        auto posJobInt = job;
        seqCop[posJobInt] = *posA;
        seqCop[IposA] = seq_swap[job];
        auto res = lbdBuildSeq(seqCop);
        if(res != std::nullopt)
            return res;

        return std::nullopt;
    }

    std::optional<SASolution> SimulatedAnnealing::move_swap_WC(const SASolution& sol)
    {
        auto mock = sol.mock;
        if (!mock) return std::nullopt;

        std::vector<unsigned short> active_ops;
        for (auto& [op_id, dec] : sol.decisions)
            if (!dec.excluded) active_ops.push_back(op_id);
        if (active_ops.size() < 2) return std::nullopt;

        std::uniform_int_distribution<int> dist(0, (int)active_ops.size() - 1);
        unsigned short op_a = active_ops[dist(randomEngine)];
        const Decision& dec_a = sol.decisions.at(op_a);

        struct SwapCandidate {
            unsigned short op_a;
            unsigned short op_b;
            unsigned int   empR_a_new; // empR que prendra op_a après le swap
            unsigned int   empR_b_new;
        };
        std::vector<SwapCandidate> swap_candidates;

        // Chercher une op dans une EmpV différent compatible
        //pair < op, empR post swap>
        for (auto& op_id : active_ops) {
            if (op_id == op_a) continue;
            const Decision& dec_b = sol.decisions.at(op_id);
            if (dec_b.empV == dec_a.empV) continue;

            // Vérifier compatibilité infrastructure
            if (!STFMockInstance::infraComp[op_a][dec_b.voie]) continue;
            if (!STFMockInstance::infraComp[op_id][dec_a.voie]) continue;
            unsigned short empRA = 0;
            unsigned short empRB = 0;
            auto findDispOpA = std::find_if(mock->jobDispoVoiesRames[op_a].begin(), mock->jobDispoVoiesRames[op_a].end(), [&](auto& el){
                return mock->dispoVoiesRames[el].dispoVoie == dec_b.empV;
            });
            auto findDispOpB = std::find_if(mock->jobDispoVoiesRames[op_id].begin(), mock->jobDispoVoiesRames[op_id].end(), [&](auto& el){
                return mock->dispoVoiesRames[el].dispoVoie == dec_a.empV;
            });
            if(findDispOpA != mock->jobDispoVoiesRames[op_a].end() && findDispOpB != mock->jobDispoVoiesRames[op_id].end())
            {
                empRA = mock->dispoVoiesRames[*findDispOpA].dispoRame;
                empRB = mock->dispoVoiesRames[*findDispOpB].dispoRame;
            }
            else
                continue;

            swap_candidates.push_back({
                op_a, op_id, empRA,empRB
            });
        }
        if (swap_candidates.empty()) return std::nullopt;

        // Tirer deuxième op parmi les candidats
        std::uniform_int_distribution<int> cand_dist(0, (int)swap_candidates.size() - 1);
        auto swap = swap_candidates[cand_dist(randomEngine)];
        const Decision& dec_b = sol.decisions.at(swap.op_b);

        // échanger tracks et slots
        SASolution neighbor = sol;
        Decision& new_dec_a = neighbor.decisions[op_a];
        Decision& new_dec_b = neighbor.decisions[swap.op_b];

        new_dec_a.voie        = dec_b.voie;
        new_dec_a.site        = dec_b.site;
        new_dec_a.empV        = dec_b.empV;
        new_dec_a.empR        = swap.empR_a_new;
        new_dec_a.timeslotGraphSplited = mock->trajectoryStops[swap.empR_a_new].getDispoStop();
        new_dec_a.lastCreneau = {0,0};

        new_dec_b.voie        = dec_a.voie;
        new_dec_b.site        = dec_a.site;
        new_dec_b.empV        = dec_a.empV;
        new_dec_b.empR        = swap.empR_b_new;
        new_dec_b.timeslotGraphSplited = mock->trajectoryStops[swap.empR_b_new].getDispoStop();
        new_dec_b.lastCreneau = {0,0};

        auto oldCrenA = dec_a.lastCreneau;
        auto oldCrenB = dec_b.lastCreneau;
        
        std::vector<std::pair<unsigned short, decision>> jobsEmpVA;
        std::vector<std::pair<unsigned short, decision>> jobsEmpVB;
        unsigned int oldCost = dec_a.lastCreneau.first * STFMockInstance::jobs[op_a]->getPoidsRetard() + dec_b.lastCreneau.first * STFMockInstance::jobs[swap.op_b]->getPoidsRetard();
        for(auto& dec : neighbor.decisions)
        {
            if(!dec.second.excluded)
            {
                if(dec.second.empV == new_dec_a.empV)
                {
                    jobsEmpVA.push_back({dec.first, dec.second});
                    oldCost += dec.second.lastCreneau.first * STFMockInstance::jobs[dec.first]->getPoidsRetard();
                }
                if(dec.second.empV == new_dec_b.empV)
                {
                    jobsEmpVB.push_back({dec.first, dec.second});
                    oldCost += dec.second.lastCreneau.first * STFMockInstance::jobs[dec.first]->getPoidsRetard();
                }
            }
        }
        std::sort(jobsEmpVA.begin(), jobsEmpVA.end(), [&](auto& el1, auto& el2){
            auto cren1 = el1.second.lastCreneau;
            auto cren2 = el2.second.lastCreneau;
            if(cren1.first == 0 && cren1.second == 0)
            {
                cren1 = oldCrenB;
            }
            if(cren2.first == 0 && cren2.second == 0)
            {
                cren2 = oldCrenB;
            }
            return cren1.first < cren2.first;
        });

        std::sort(jobsEmpVB.begin(), jobsEmpVB.end(), [&](auto& el1, auto& el2){
            auto cren1 = el1.second.lastCreneau;
            auto cren2 = el2.second.lastCreneau;
            if(cren1.first == 0 && cren1.second == 0)
            {
                cren1 = oldCrenA;
            }
            if(cren2.first == 0 && cren2.second == 0)
            {
                cren2 = oldCrenA;
            }
            return cren1.first < cren2.first;
        });
        auto resA = checkSequence(jobsEmpVA, new_dec_a.empV);
        auto resB = checkSequence(jobsEmpVB, new_dec_b.empV);
        unsigned int newCost = 0;
        if(resA && resB)
        {
            unsigned int id = 0;
            for(auto& jobsA : jobsEmpVA)
            {
                auto& dec = neighbor.decisions[jobsA.first];
                dec.lastCreneau = resA.value()[id].second.lastCreneau;
                newCost += dec.lastCreneau.first * STFMockInstance::jobs[jobsA.first]->getPoidsRetard();
                ++id;
            }
            id = 0;
            for(auto& jobsB : jobsEmpVB)
            {
                auto& dec = neighbor.decisions[jobsB.first];
                dec.lastCreneau = resB.value()[id].second.lastCreneau;
                newCost += dec.lastCreneau.first * STFMockInstance::jobs[jobsB.first]->getPoidsRetard();
                ++id;
            }
        }
        else
            return std::nullopt;
        neighbor.cost = (neighbor.cost - oldCost) + newCost;
        neighbor.source = ESourceTrackPlan::SimAn;
        return neighbor;
    }

    std::optional<SASolution> SimulatedAnnealing::move_insert_WC(const SASolution& sol)
    {
        auto mock = sol.mock;
        if (!mock) return std::nullopt;

        std::vector<unsigned short> inactive_ops;
        for (auto& [op_id, dec] : sol.decisions)
            if (dec.excluded) inactive_ops.push_back(op_id);
        if (inactive_ops.empty()) return std::nullopt;

        std::uniform_int_distribution<int> dist(0, (int)inactive_ops.size() - 1);
        unsigned short op_a = inactive_ops[dist(randomEngine)];

        std::vector<unsigned int> dispCandidate;
        for (auto& disp : mock->jobDispoVoiesRames[op_a]) {

            auto dur = mock->dispoVoiesRames[disp].match.second - mock->dispoVoiesRames[disp].match.first;
            if(dur >= STFMockInstance::jobs[op_a]->getDuree())
            {
                dispCandidate.push_back(disp);
            }
            else if(dur >= STFMockInstance::jobs[op_a]->getDureeDiag() && sol.diagCost + STFMockInstance::jobs[op_a]->getPoidsRejet() <= configlib::Configuration::Global.EPSILON)
            {
                dispCandidate.push_back(disp);
            }
        }
        if(dispCandidate.empty()) return std::nullopt;

        std::uniform_int_distribution<int> cand_dist(0, (int)dispCandidate.size() - 1);
        auto disp = dispCandidate[cand_dist(randomEngine)];

        //Construire le voisin - try insert
        SASolution neighbor = sol;
        Decision& new_dec_a = neighbor.decisions[op_a];
        auto dur = mock->dispoVoiesRames[disp].match.second - mock->dispoVoiesRames[disp].match.first;

        new_dec_a.rejected = dur >= STFMockInstance::jobs[op_a]->getDuree() ? false : true,
        new_dec_a.excluded = false;
        new_dec_a.voie        = mock->dispoVoiesRames[disp].voie;
        new_dec_a.site        = mock->dispoVoiesRames[disp].site;
        new_dec_a.empV        = mock->dispoVoiesRames[disp].dispoVoie;
        new_dec_a.empR        = mock->dispoVoiesRames[disp].dispoRame;
        new_dec_a.timeslotGraphSplited = mock->trajectoryStops[new_dec_a.empR].getDispoStop();
        new_dec_a.lastCreneau = {0,0};
        
        unsigned int oldCost = MAXIMUM_TIME_OFFSET*STFMockInstance::jobs[op_a]->getPoidsRetard();
        unsigned int oldDiagCost = 0;
        unsigned int newDiagCost = new_dec_a.rejected ? STFMockInstance::jobs[op_a]->getPoidsRejet() : 0;
        std::vector<std::pair<unsigned short, decision>> jobsEmpVA;
        for(auto& dec : neighbor.decisions)
        {
            if(!dec.second.excluded && dec.first != op_a)
            {
                if(dec.second.empV == new_dec_a.empV)
                {
                    jobsEmpVA.push_back({dec.first, dec.second});
                    oldCost += dec.second.lastCreneau.first * STFMockInstance::jobs[dec.first]->getPoidsRetard();
                }
            }
        }

        std::sort(jobsEmpVA.begin(), jobsEmpVA.end(), [&](auto& el1, auto& el2){
            auto cren1 = el1.second.lastCreneau;
            auto cren2 = el2.second.lastCreneau;
            return cren1.first < cren2.first;
        });

        std::uniform_int_distribution<unsigned int> distPos(0,jobsEmpVA.size());
        auto pos = distPos(randomEngine);
        
        auto seq = jobsEmpVA;
        if(pos == jobsEmpVA.size())
        {
            seq.insert(seq.end(), std::make_pair(op_a, new_dec_a));
        }
        else {
            seq.insert(seq.begin() + pos, std::make_pair(op_a, new_dec_a));
        }
        unsigned int newCost = 0;
        auto res = checkSequence(seq, new_dec_a.empV);
        if(res)
        {
            unsigned int id = 0;
            for(auto& jobsA : seq)
            {
                auto& dec = neighbor.decisions[jobsA.first];
                dec.lastCreneau = res.value()[id].second.lastCreneau;
                newCost += dec.lastCreneau.first * STFMockInstance::jobs[jobsA.first]->getPoidsRetard();
                ++id;
            }
            neighbor.diagCost = (neighbor.diagCost - oldDiagCost) + newDiagCost;
            neighbor.cost = (neighbor.cost - oldCost) + newCost;
            neighbor.source = ESourceTrackPlan::SimAn;
            return neighbor;
        }

        return std::nullopt;
    }

    std::optional<SASolution> SimulatedAnnealing::move_move_WC(const SASolution& sol)
    {
        auto mock = sol.mock;
        if (!mock) return std::nullopt;

        std::vector<unsigned short> active_ops;
        for (auto& [op_id, dec] : sol.decisions)
            if (!dec.excluded) active_ops.push_back(op_id);
        if (active_ops.empty()) return std::nullopt;

        std::uniform_int_distribution<int> dist(0, (int)active_ops.size() - 1);
        unsigned short op_a = active_ops[dist(randomEngine)];

        std::vector<unsigned int> dispCandidate;
        const Decision& dec_a = sol.decisions.at(op_a);

        for (auto& disp : mock->jobDispoVoiesRames[op_a]) {
            if(dec_a.empV == mock->dispoVoiesRames[disp].dispoVoie)
                continue;

            auto dur = mock->dispoVoiesRames[disp].match.second - mock->dispoVoiesRames[disp].match.first;
            if(dur >= STFMockInstance::jobs[op_a]->getDuree())
            {
                dispCandidate.push_back(disp);
            }
            else if(dur >= STFMockInstance::jobs[op_a]->getDureeDiag() && sol.diagCost + STFMockInstance::jobs[op_a]->getPoidsRejet() <= configlib::Configuration::Global.EPSILON)
            {
                dispCandidate.push_back(disp);
            }
        }
        if(dispCandidate.empty()) return std::nullopt;

        std::uniform_int_distribution<int> cand_dist(0, (int)dispCandidate.size() - 1);
        auto disp = dispCandidate[cand_dist(randomEngine)];

        //Construire le voisin - try insert
        SASolution neighbor = sol;
        Decision& new_dec_a = neighbor.decisions[op_a];
        auto dur = mock->dispoVoiesRames[disp].match.second - mock->dispoVoiesRames[disp].match.first;

        new_dec_a.rejected = dur >= STFMockInstance::jobs[op_a]->getDuree() ? false : true,
        new_dec_a.voie        = mock->dispoVoiesRames[disp].voie;
        new_dec_a.site        = mock->dispoVoiesRames[disp].site;
        new_dec_a.empV        = mock->dispoVoiesRames[disp].dispoVoie;
        new_dec_a.empR        = mock->dispoVoiesRames[disp].dispoRame;
        new_dec_a.timeslotGraphSplited = mock->trajectoryStops[new_dec_a.empR].getDispoStop();
        new_dec_a.lastCreneau = {0,0};
        
        unsigned int oldCost = dec_a.lastCreneau.first*STFMockInstance::jobs[op_a]->getPoidsRetard();
        unsigned int oldDiagCost = dec_a.rejected ? STFMockInstance::jobs[op_a]->getPoidsRejet() : 0;
        unsigned int newDiagCost = new_dec_a.rejected ? STFMockInstance::jobs[op_a]->getPoidsRejet() : 0;

        std::vector<std::pair<unsigned short, decision>> jobsEmpVA;
        std::vector<std::pair<unsigned short, decision>> jobsEmpVB;

        for(auto& dec : neighbor.decisions)
        {
            if(!dec.second.excluded && dec.first != op_a)
            {
                if(dec.second.empV == new_dec_a.empV)
                {
                    jobsEmpVA.push_back({dec.first, dec.second});
                    oldCost += dec.second.lastCreneau.first * STFMockInstance::jobs[dec.first]->getPoidsRetard();
                }
            }
            if(!dec.second.excluded)
            {
                if(dec.second.empV == dec_a.empV)
                {
                    jobsEmpVB.push_back({dec.first, dec.second});
                    oldCost += dec.second.lastCreneau.first * STFMockInstance::jobs[dec.first]->getPoidsRetard();
                }
            }
        }

        std::sort(jobsEmpVA.begin(), jobsEmpVA.end(), [&](auto& el1, auto& el2){
            auto cren1 = el1.second.lastCreneau;
            auto cren2 = el2.second.lastCreneau;
            return cren1.first < cren2.first;
        });
        std::sort(jobsEmpVB.begin(), jobsEmpVB.end(), [&](auto& el1, auto& el2){
            auto cren1 = el1.second.lastCreneau;
            auto cren2 = el2.second.lastCreneau;
            return cren1.first < cren2.first;
        });

        unsigned int newCost = 0;
        //décale à gauche sur track de départ
        auto resB = checkSequence(jobsEmpVB, dec_a.empV);
        unsigned int id = 0;
        for(auto& jobsB : jobsEmpVB)
        {
            auto& dec = neighbor.decisions[jobsB.first];
            dec.lastCreneau = resB.value()[id].second.lastCreneau;
            newCost+= dec.lastCreneau.first * STFMockInstance::jobs[jobsB.first]->getPoidsRetard();
            ++id;
        }
        

        std::uniform_int_distribution<unsigned int> distPos(0,jobsEmpVA.size());
        auto pos = distPos(randomEngine);
        
        auto seq = jobsEmpVA;
        if(pos == jobsEmpVA.size())
        {
            seq.insert(seq.end(), std::make_pair(op_a, new_dec_a));
        }
        else {
            seq.insert(seq.begin() + pos, std::make_pair(op_a, new_dec_a));
        }
        auto res = checkSequence(seq, new_dec_a.empV);
        if(res)
        {
            unsigned int id = 0;
            for(auto& jobsA : seq)
            {
                auto& dec = neighbor.decisions[jobsA.first];
                dec.lastCreneau = res.value()[id].second.lastCreneau;
                newCost+= dec.lastCreneau.first * STFMockInstance::jobs[jobsA.first]->getPoidsRetard();
                ++id;
            }
            neighbor.diagCost = (neighbor.diagCost - oldDiagCost) + newDiagCost;
            neighbor.cost = (neighbor.cost - oldCost) + newCost;
            neighbor.source = ESourceTrackPlan::SimAn;
            return neighbor;
        }

        return std::nullopt;
    }

    std::optional<SASolution> SimulatedAnnealing::move_remove_WC(const SASolution& sol)
    {
        auto mock = sol.mock;
        if (!mock) return std::nullopt;

        std::vector<unsigned short> active_ops;
        for (auto& [op_id, dec] : sol.decisions)
            if (!dec.excluded) active_ops.push_back(op_id);
        if (active_ops.empty()) return std::nullopt;

        std::uniform_int_distribution<int> dist(0, (int)active_ops.size() - 1);
        unsigned short op_a = active_ops[dist(randomEngine)];
        //const Decision& dec_a = sol.decisions.at(op_a);

        //Construire le voisin - exclure a
        SASolution neighbor = sol;
        const Decision& dec_a = sol.decisions.at(op_a);
        Decision& new_dec_a = neighbor.decisions[op_a];

        new_dec_a.rejected = false,
        new_dec_a.excluded = true;
        new_dec_a.voie        = 0;
        new_dec_a.site        = 0;
        new_dec_a.empV        = 0;
        new_dec_a.empR        = 0;
        new_dec_a.timeslotGraphSplited = CreneauHoraire();
        new_dec_a.lastCreneau = {0,0};
        
        unsigned int oldCost = dec_a.lastCreneau.first * STFMockInstance::jobs[op_a]->getPoidsRetard();
        unsigned int oldDiagCost = dec_a.rejected ? STFMockInstance::jobs[op_a]->getPoidsRejet() : 0;
        unsigned int newDiagCost = 0;

        std::vector<std::pair<unsigned short, decision>> jobsEmpVA;
        for(auto& dec : neighbor.decisions)
        {
            if(!dec.second.excluded)
            {
                if(dec.second.empV == dec_a.empV)
                {
                    jobsEmpVA.push_back({dec.first, dec.second});
                    oldCost += dec.second.lastCreneau.first * STFMockInstance::jobs[dec.first]->getPoidsRetard();
                }
            }
        }
        unsigned int newCost = MAXIMUM_TIME_OFFSET * STFMockInstance::jobs[op_a]->getPoidsRetard();

        if(!jobsEmpVA.empty())
        {

            std::sort(jobsEmpVA.begin(), jobsEmpVA.end(), [&](auto& el1, auto& el2){
                auto cren1 = el1.second.lastCreneau;
                auto cren2 = el2.second.lastCreneau;
                return cren1.first < cren2.first;
            });
            auto resA = checkSequence(jobsEmpVA, dec_a.empV);
            if(resA)
            {
                unsigned int id = 0;
                for(auto& jobsA : jobsEmpVA)
                {
                    auto& dec = neighbor.decisions[jobsA.first];
                    dec.lastCreneau = resA.value()[id].second.lastCreneau;
                    newCost += dec.lastCreneau.first * STFMockInstance::jobs[jobsA.first]->getPoidsRetard();
                    ++id;
                }
            }
            else
                return std::nullopt;
        }
        neighbor.diagCost = (neighbor.diagCost - oldDiagCost) + newDiagCost;
        neighbor.cost = (neighbor.cost - oldCost) + newCost;
        neighbor.source = ESourceTrackPlan::SimAn;
        return neighbor;
    }

    std::optional<SASolution> SimulatedAnnealing::move_change_mode_WC(const SASolution& sol)
    {
        auto mock = sol.mock;
        if (!mock) return std::nullopt;

        std::vector<unsigned short> active_ops;
        for (auto& [op_id, dec] : sol.decisions)
            if (!dec.excluded) active_ops.push_back(op_id);
        if (active_ops.empty()) return std::nullopt;

        std::uniform_int_distribution<int> dist(0, (int)active_ops.size() - 1);
        unsigned short op_a = active_ops[dist(randomEngine)];
        const Decision& dec_a = sol.decisions.at(op_a);

        SASolution neighbor = sol;
        Decision& new_dec_a = neighbor.decisions[op_a];

        unsigned int oldDiagCost = dec_a.rejected ? STFMockInstance::jobs[op_a]->getPoidsRejet() : 0;
        unsigned int newDiagCost = dec_a.rejected ? 0 : STFMockInstance::jobs[op_a]->getPoidsRejet();
        unsigned int oldCost = 0;
        if(dec_a.rejected)
        {
            new_dec_a.rejected = false;
        }
        else {
            if(neighbor.diagCost + STFMockInstance::jobs[op_a]->getPoidsRejet() <= configlib::Configuration::Global.EPSILON)
            {
                new_dec_a.rejected = true;
            }
            else {
                return std::nullopt;
            }
        }
        
        std::vector<std::pair<unsigned short, decision>> jobsEmpVA;
        for(auto& dec : neighbor.decisions)
        {
            if(!dec.second.excluded)
            {
                if(dec.second.empV == dec_a.empV)
                {
                    jobsEmpVA.push_back({dec.first, dec.second});
                    oldCost += dec.second.lastCreneau.first * STFMockInstance::jobs[dec.first]->getPoidsRetard();
                }
            }
        }
        
        std::sort(jobsEmpVA.begin(), jobsEmpVA.end(), [&](auto& el1, auto& el2){
            auto cren1 = el1.second.lastCreneau;
            auto cren2 = el2.second.lastCreneau;
            return cren1.first < cren2.first;
        });
        unsigned int newCost = 0;
        if(!jobsEmpVA.empty())
        {
            auto resA = checkSequence(jobsEmpVA, dec_a.empV);
            if(resA)
            {
                unsigned int id = 0;
                for(auto& jobsA : jobsEmpVA)
                {
                    auto& dec = neighbor.decisions[jobsA.first];
                    dec.lastCreneau = resA.value()[id].second.lastCreneau;
                    newCost += dec.lastCreneau.first * STFMockInstance::jobs[jobsA.first]->getPoidsRetard();
                    ++id;
                }
            }
            else
                return std::nullopt;
        }
        neighbor.diagCost = (neighbor.diagCost - oldDiagCost) + newDiagCost;
        neighbor.cost = (neighbor.cost - oldCost) + newCost;
        neighbor.source = ESourceTrackPlan::SimAn;
        return neighbor;
    }

    std::optional<std::vector<std::pair<unsigned short, Decision>>> SimulatedAnnealing::checkSequence(const std::vector<std::pair<unsigned short, Decision>>& jobsDec, unsigned int machine)
    {
        auto machineDisp = STFMockInstance::machines[machine]->getDispo();
        unsigned short minBegin = machineDisp.getDebut().getRelativeDate();
        std::vector<std::pair<unsigned short, Decision>> res;
        for(auto& jobDec : jobsDec)
        {
            auto matchWithMachine = CreneauHoraire::checkSlotsCompatibility(jobDec.second.timeslotGraphSplited, machineDisp).second;
            minBegin = std::max(minBegin, matchWithMachine.first);
            auto newJobDec = jobDec;
            unsigned int duration = !jobDec.second.rejected ? STFMockInstance::jobs[jobDec.first]->getDuree() : STFMockInstance::jobs[jobDec.first]->getDureeDiag();
            newJobDec.second.lastCreneau = {minBegin, minBegin + duration};
            if(minBegin + duration > matchWithMachine.first + matchWithMachine.second)
                return std::nullopt;
            minBegin = minBegin + duration;
            res.push_back(newJobDec);
        }
        return std::optional<std::vector<std::pair<unsigned short, Decision>>>(res);
    }


    std::pair<bool, std::vector<unsigned int>> SimulatedAnnealing::PSE_Carlier_Rivreau(std::vector<std::pair<unsigned short, decision>> &jobs, unsigned int voieMachine) {

        std::stringstream fakeFile;

        for (unsigned int i = 0; i < jobs.size(); ++i) {
            auto match = CreneauHoraire::checkSlotsCompatibility(jobs[i].second.timeslotGraphSplited,
                                                                    STFMockInstance::machines[voieMachine]->getDispo());
            fakeFile << i + 1;
            fakeFile << " " << match.second.first;
            if (!jobs[i].second.rejected) {
                fakeFile << " " << STFMockInstance::jobs[jobs[i].first]->getDuree();
            } else
                fakeFile << " " << STFMockInstance::jobs[jobs[i].first]->getDureeDiag();
            fakeFile << " " << match.second.first + match.second.second;
            fakeFile << std::endl;
        }
        std::ifstream file;
        file.basic_ios<char>::rdbuf(fakeFile.rdbuf());
        Solution sol((int) jobs.size());
        OneMachine machine(file);
        machine.solve(sol);

        std::vector<unsigned int> solution(jobs.size());
        for (unsigned int i = 0; i < jobs.size(); ++i) {
            solution[i] = sol.startTime[i+1];
        }

        return {machine.checkSol(sol) && sol.Lmax <= 0, solution};
    }


    /*unsigned long SimulatedAnnealing::getUniquePlans(std::vector<TrackPlan>& plans)
    {
        std::set<TrackPlan> uniqueSchedules;
        std::vector<TrackPlan> newTMPVec;
        for(auto& trackSch : plans)
        {
            uniqueSchedules.insert(trackSch);
        }
        auto uniqueNb = uniqueSchedules.size();
        std::move(uniqueSchedules.begin(), uniqueSchedules.end(), std::back_inserter(newTMPVec));
        std::swap(plans, newTMPVec);

        return uniqueNb;
    }*/



}