CMakeLists.txt

@@ -28,6 +28,7 @@ add_executable(pdptw src/mains/main.cpp
                 src/lns/solution/solution.cpp
                 src/lns/constraints/constraint.cpp
                 src/lns/constraints/capacity/capacity_constraint.cpp
+                src/lns/constraints/time_window/forward_time_slack.cpp
                 src/lns/constraints/time_window/time_window_constraint.cpp
                 src/lns/modification/pair/insert_pair.cpp 
                 src/lns/modification/pair/remove_pair.cpp
@@ -38,7 +39,9 @@ add_executable(pdptw src/mains/main.cpp
                 src/output/solution_checker.cpp
                 src/output/solution_exporter.cpp
                 src/lns/operators/sorting_strategy.cpp
+                src/lns/operators/destruction/clean_empty_route.cpp
                 src/lns/operators/destruction/random_destroy.cpp
+                src/lns/operators/destruction/string_removal.cpp
                 src/lns/operators/reconstruction/enumerate.cpp
                 src/lns/operators/reconstruction/list_heuristic_insertion.cpp
                 src/lns/operators/reconstruction/list_heuristic_cost_oriented.cpp

src/config.h

-
-
 const int EXCLUSION_PENALTY = 100;
 
-const int RANDOM_SEED = 100;
+const int RANDOM_SEED = 100;
\ No newline at end of file
+
+const int NUMBER_VEHICLE = 20;
\ No newline at end of file

src/input/data.cpp

@@ -24,23 +24,23 @@ double data::routeCost(PDPTWData const & data, Route const & route)
 {
     const Matrix & matrix = data.getMatrix();
     const std::vector<int> & routeIDs = route.getRoute();
-    double cost = 0;
+    double cost = 0.0;
 
     if (routeIDs.empty())
     {
-        return 0;
+        return 0.0;
     }
     // cost from and to the depot
     cost +=  matrix.at(0).at(routeIDs.at(0));
-    //std::cout << "\n route cost : " << matrix.at(0).at(routeIDs.at(0) << " ";
+    // std::cout << "\n route cost : " << matrix.at(0).at(routeIDs.at(0)) << " ";
     cost +=  matrix.at(routeIDs.back()).at(0);
     
     // cost in the route
     for (size_t i = 0; i < routeIDs.size() - 1; ++i) {
         cost += matrix.at(routeIDs.at(i)).at(routeIDs.at(i+1));
-        //std::cout << matrix.at(routeIDs.at(i).at(routeIDs.at(i+1) << " ";
+        // std::cout << matrix.at(routeIDs.at(i)).at(routeIDs.at(i+1)) << " ";
     }
-    //std::cout << matrix.at(routeIDs.back()).at(0) << " : " << cost << "\n";
+    // std::cout << matrix.at(routeIDs.back()).at(0) << " : " << cost << "\n";
 
     return cost;
 }

src/input/pdptw_data.cpp

 #include "pdptw_data.h"
 
-#include <iostream>
+#include "data.h"
+
+#include <algorithm>
 #include <fstream>
+#include <iostream>
+#include <numeric>
 #include <spdlog/spdlog.h>
 
-
 int PDPTWData::getSize() const
 {
     return size;
@@ -15,12 +18,12 @@ int PDPTWData::getCapacity() const
     return capacity;
 }
 
-std::vector<Location> const &PDPTWData::getLocations() const 
+std::vector<Location> const &PDPTWData::getLocations() const
 {
     return locations;
 }
 
-std::vector<Pair> const &PDPTWData::getPairs() const 
+std::vector<Pair> const &PDPTWData::getPairs() const
 {
     return pairs;
 }
@@ -37,12 +40,12 @@ std::string PDPTWData::getDataName() const
 
 Location const &PDPTWData::getLocation(int id) const
 {
-    if (id==0)
+    if (id == 0)
     {
         return getDepot();
     }
     // location index from 0 to n-1
-    return locations.at(id -1);
+    return locations.at(id - 1);
 }
 
 Matrix const &PDPTWData::getMatrix() const
@@ -50,24 +53,28 @@ Matrix const &PDPTWData::getMatrix() const
     return distanceMatrix;
 }
 
-PDPTWData::PDPTWData(std::string dataName, int size, int capacity, Location depot, std::vector<Location> locations, Matrix distanceMatrix)
-    : dataName(dataName), size(size), capacity(capacity), depot(depot), locations(std::move(locations)), distanceMatrix(std::move(distanceMatrix)) 
+PDPTWData::PDPTWData(std::string dataName, int size, int capacity, Location depot, std::vector<Location> locations,
+                     Matrix distanceMatrix)
+    : dataName(dataName), size(size), capacity(capacity), depot(depot), locations(std::move(locations)),
+      distanceMatrix(std::move(distanceMatrix))
 {
     // Associate pair of locations
     pairs.clear();
-    for (const Location & loc : this->locations)
+    for (Location const &loc: this->locations)
     {
-        if( loc.getLocType() == LocType::PICKUP )
+        if (loc.getLocType() == LocType::PICKUP)
         {
             // vector indexed from 0 / Location indexed from 1
-            pairs.emplace_back(loc, this->locations.at(loc.getPair()-1), loc.getId());
+            pairs.emplace_back(loc, this->locations.at(loc.getPair() - 1), loc.getId());
         }
     }
+    // Compute closest location matrix
+    initClosestLocations();
 }
 
-const Pair &PDPTWData::getPair(int id) const
+Pair const &PDPTWData::getPair(int id) const
 {
-    for (const Pair &pair : pairs)
+    for (Pair const &pair: pairs)
     {
         if (id == pair.getID())
         {
@@ -83,8 +90,12 @@ int PDPTWData::getPairCount() const
     return getPairs().size();
 }
 
+int PDPTWData::getLocationCount() const
+{
+    return getLocations().size();
+}
 
-void PDPTWData::print() const 
+void PDPTWData::print() const
 {
     std::cout << "Instance name : " << dataName << "\n";
     std::cout << "Instance size: " << size << "\n";
@@ -93,22 +104,36 @@ void PDPTWData::print() const
     depot.print();
 
     std::cout << "Locations:\n";
-    for (const auto& loc : locations) {
+    for (auto const &loc: locations)
+    {
         loc.print();
     }
 
     std::cout << "Distance Matrix:\n";
-    for (const auto& row : distanceMatrix) {
-        for (const auto& dist : row) {
+    for (auto const &row: distanceMatrix)
+    {
+        for (auto const &dist: row)
+        {
             std::cout << dist << " ";
         }
         std::cout << "\n";
     }
 
+    std::cout << "Closest Matrix (no depot):\n";
+    for (auto const &row: closestLocations)
+    {
+        for (auto const &ID: row)
+        {
+            std::cout << ID << " ";
+        }
+        std::cout << "\n";
+    }
+
     std::cout << "Pair IDs:\n";
-    for (const auto& pair : pairs)
-    {;
-        std::cout << pair.getID() << " "; 
+    for (auto const &pair: pairs)
+    {
+        ;
+        std::cout << pair.getID() << " ";
     }
     std::cout << " \n";
 }
@@ -127,14 +152,13 @@ void PDPTWData::checkData() const
     }
 }
 
-
 bool PDPTWData::checkMatrix() const
 {
     // square matrix
-    
-    for (const auto& row : getMatrix()) 
+
+    for (auto const &row: getMatrix())
     {
-        if (row.size() != size) 
+        if (row.size() != size)
         {
             return true;
         }
@@ -161,22 +185,23 @@ bool PDPTWData::checkMatrix() const
 
 bool PDPTWData::checkLocation() const
 {
-
     // check if location id equals the position in the location vector
-    for (size_t i = 0; i < size-1; ++i) {
-        if (locations.at(i).getId() != static_cast<int>(i)+1) {
+    for (size_t i = 0; i < size - 1; ++i)
+    {
+        if (locations.at(i).getId() != static_cast<int>(i) + 1)
+        {
             return true;
         }
     }
-    
+
     // check if pair of location are well made (type, id, demand, timeWindow)
-    for(const Location& loc : getLocations()) 
+    for (Location const &loc: getLocations())
     {
         if (loc.getLocType() == LocType::PICKUP)
-        {   
-            if ( (getLocations().at(loc.getPair()-1).getLocType() != LocType::DELIVERY) 
-            || (loc.getDemand() != - getLocations().at(loc.getPair()-1).getDemand()) 
-            || (loc.getId() != getLocations().at(loc.getPair()-1).getPair()) )
+        {
+            if ((getLocations().at(loc.getPair() - 1).getLocType() != LocType::DELIVERY) ||
+                (loc.getDemand() != -getLocations().at(loc.getPair() - 1).getDemand()) ||
+                (loc.getId() != getLocations().at(loc.getPair() - 1).getPair()))
             {
                 return true;
             }
@@ -191,6 +216,40 @@ bool PDPTWData::checkLocation() const
     return false;
 }
 
+void PDPTWData::initClosestLocations()
+{
+    closestLocations.reserve(getLocationCount());
+    for (Location const &location: locations)
+    {
+        // sorting from closest to furthest using indexes
+        std::vector<int> closestLocationsIndexes(getLocationCount());
+        std::iota(closestLocationsIndexes.begin(), closestLocationsIndexes.end(), 0);
+
+        // Do the sorting
+        std::ranges::sort(closestLocationsIndexes, {}, [this, &location](int index) {
+            return data::TravelTime(*this, location.getId(), locations.at(index).getId());
+        });
+
+        // we store the indexes in the final container
+        closestLocations.emplace_back();
+        closestLocations.back().reserve(getLocationCount());
+        std::ranges::transform(closestLocationsIndexes, std::back_inserter(closestLocations.back()), [this](int index) {
+            return std::cref(locations.at(index)).get().getId();
+        });
+    }
+}
+
+std::vector<int> const &PDPTWData::getClosestLocationsID(int id) const
+{
+    // sanity check
+    if (id < 1 || id > getLocationCount())
+    {
+        spdlog::error("getClosestLocationsID: Location ID {} is out of valid range [1, {}]", id, getLocationCount());
+        throw std::out_of_range("Error: Location ID is out of valid range [1, " + std::to_string(getLocationCount()) +
+                                "]");
+    }
+    return closestLocations.at(id - 1);// vector indexed from 0
+}
 
 InputJsonException::InputJsonException(std::string_view reason)
     : reason(fmt::format("Input JSON file is incorrect : {}", reason))

src/input/pdptw_data.h

@@ -30,6 +30,11 @@ class PDPTWData
     std::vector<Location> locations;
     std::vector<Pair> pairs;// std::unordered_map<int, Pair> pair; if getPair(index) is needed ?
     Matrix distanceMatrix;
+    // stores the list of locations ID from closest to furthest for every location (no depot).
+    std::vector<std::vector<int>> closestLocations;
+
+    void initClosestLocations();
+
 
 public:
     PDPTWData();
@@ -54,19 +59,25 @@ public:
 
 
     std::vector<Location> const &getLocations() const;
+    int getLocationCount() const;
     std::vector<Pair> const &getPairs() const;
     int getPairCount() const;
     Pair const &getPair(int id) const;
 
     /** 
-    *   0 return the depot.
-    *   Other numbers return the associated location.
-    */
+     *  0 return the depot.
+     *  Other numbers return the associated location.
+     */
     Location const &getLocation(int id) const;
     Location const &getDepot() const;
 
     Matrix const &getMatrix() const;
 
+    /** 
+     *  Given a location ID, return a list of location id sorted by proximity.
+     */
+    std::vector<int> const &getClosestLocationsID(int id) const;
+
     int getSize() const;
     int getCapacity() const;
     std::string getDataName() const;

src/lns/constraints/capacity/capacity_constraint.cpp

 #include "capacity_constraint.h"
 
-
 CapacityConstraint::CapacityConstraint(Solution const &solution) : Constraint(solution)
 {
     // Init an empty vector of reach time for all routes
-    //std::vector<CapacityVector> routeCapacities(getSolution().getRoutes().size(), CapacityVector());
     routeCapacities.resize(getSolution().getRoutes().size());
 }
 
-std::unique_ptr<Constraint> CapacityConstraint::clone(Solution const &newOwningSolution) const {
+std::unique_ptr<Constraint> CapacityConstraint::clone(Solution const &newOwningSolution) const
+{
     std::unique_ptr<CapacityConstraint> clonePtr = std::make_unique<CapacityConstraint>(newOwningSolution);
     clonePtr->routeCapacities = routeCapacities;
     return clonePtr;
 }
 
-
-void CapacityConstraint::initCapacities()
+void CapacityConstraint::initRouteCapacities()
 {
-    routeCapacities = std::vector<CapacityVector>();
-    for (const Route& route : getSolution().getRoutes()) 
+    routeCapacities.clear();
+    for (Route const &route: getSolution().getRoutes())
     {
         CapacityVector capacities;
         int currentCapacity = 0;
-        
-        for (int locationID : route.getRoute()) 
+
+        for (int locationID: route.getRoute())
         {
             currentCapacity += getSolution().getData().getLocation(locationID).getDemand();
             capacities.push_back(currentCapacity);
         }
-        
+
         routeCapacities.push_back(capacities);
     }
 }
 
 
-std::vector<int> const & CapacityConstraint::getRouteCapacities(int routeIndex) const
+std::vector<int> const &CapacityConstraint::getRouteCapacities(int routeIndex) const
 {
     return routeCapacities.at(routeIndex);
 }
 
-// check for every location between the pickupPosition and the deliveryPosition
-// not ideal
+/*
+ *  Check if an insertion is feasible by computing new capacity values between the PickupPosition and DeliveryPosition.
+ *  Assume the constraint is valid beforehand.
+ *  O(n), n the size of the route 
+ */
 bool CapacityConstraint::checkModif(Pair const &pair, int routeIndex, int PickupPosition, int DeliveryPosition) const
 {
     int max_capacity = getSolution().getData().getCapacity();
+    int demand = pair.getPickup().getDemand();
+    CapacityVector const &routeCapacity = getRouteCapacities(routeIndex);
 
-    // guardrail
-    if (pair.getPickup().getDemand() > max_capacity)
+    // if PickupPosition != DeliveryPosition, check from pickup-1 to delivery-1
+    // if PickupPosition is 0, don't check -1 ...
+    for (int i = std::max(0, PickupPosition - 1); i < DeliveryPosition; i++)
     {
-        return false;
-    }
-
-    // check if the disponible capacity is enough are not
-    for (int i = PickupPosition; i <= DeliveryPosition; ++i) {
-        
-        if ((i!=0) && (getRouteCapacities(routeIndex).at(i-1) + pair.getPickup().getDemand() >= max_capacity))
+        if (routeCapacity.at(i) + demand > max_capacity)
         {
             return false;
         }
-    } 
+    }
     return true;
 }
 
 // update the route and the weight
-// not ideal
-void CapacityConstraint::applyModif(Pair const &pair, int routeIndex, int PickupPosition,  int DeliveryPosition, bool addPair)
+void CapacityConstraint::applyModif(Pair const &pair, int routeIndex, int PickupPosition, int DeliveryPosition,
+                                    bool addPair)
 {
     if (addPair)
-    { 
+    {
         // Insert new values
-        routeCapacities.at(routeIndex).insert(routeCapacities.at(routeIndex).begin()+DeliveryPosition, 0);
+        routeCapacities.at(routeIndex).insert(routeCapacities.at(routeIndex).begin() + DeliveryPosition, 0);
         if (DeliveryPosition != 0)
         {
-            routeCapacities.at(routeIndex).at(DeliveryPosition) += routeCapacities.at(routeIndex).at(DeliveryPosition-1);
+            routeCapacities.at(routeIndex).at(DeliveryPosition) +=
+                    routeCapacities.at(routeIndex).at(DeliveryPosition - 1);
         }
-        routeCapacities.at(routeIndex).insert(routeCapacities.at(routeIndex).begin()+PickupPosition, pair.getPickup().getDemand());
+        routeCapacities.at(routeIndex)
+                .insert(routeCapacities.at(routeIndex).begin() + PickupPosition, pair.getPickup().getDemand());
         if (PickupPosition != 0)
         {
-            routeCapacities.at(routeIndex).at(PickupPosition) += routeCapacities.at(routeIndex).at(PickupPosition-1);
+            routeCapacities.at(routeIndex).at(PickupPosition) += routeCapacities.at(routeIndex).at(PickupPosition - 1);
         }
 
         // Update value
-        for (int i = PickupPosition + 1; i < DeliveryPosition + 1; ++i) 
+        for (int i = PickupPosition + 1; i < DeliveryPosition + 1; ++i)
         {
             routeCapacities.at(routeIndex).at(i) += pair.getPickup().getDemand();
-        } 
+        }
     }
-    else 
+    else
     {
-        for (int i = PickupPosition + 1; i < DeliveryPosition; ++i) 
+        for (int i = PickupPosition + 1; i < DeliveryPosition; ++i)
         {
             routeCapacities.at(routeIndex).at(i) += pair.getDelivery().getDemand();
         }
@@ -96,48 +96,48 @@ void CapacityConstraint::applyModif(Pair const &pair, int routeIndex, int Pickup
         // remove pair
         routeCapacities.at(routeIndex).erase(routeCapacities.at(routeIndex).begin() + DeliveryPosition);
         routeCapacities.at(routeIndex).erase(routeCapacities.at(routeIndex).begin() + PickupPosition);
-    }    
+    }
 }
 
-
 bool CapacityConstraint::check(InsertPair const &op) const
 {
     //std::cout << " #Capa Check";
     return checkModif(op.getPair(), op.getRouteIndex(), op.getPickupInsertion(), op.getDeliveryInsertion());
 }
+
 void CapacityConstraint::apply(InsertPair const &op)
 {
     //std::cout << "-> Apply Modification on Capacity \n";
     applyModif(op.getPair(), op.getRouteIndex(), op.getPickupInsertion(), op.getDeliveryInsertion(), true);
 }
 
-
 bool CapacityConstraint::check(InsertRoute const &op) const
 {
     return true;
 }
+
 void CapacityConstraint::apply(InsertRoute const &op)
 {
     // add new empty route
     routeCapacities.emplace_back();
-}   
-
+}
 
 bool CapacityConstraint::check(RemovePair const &op) const
 {
     // Remove a pair (always true)
     return true;
 }
+
 void CapacityConstraint::apply(RemovePair const &op)
 {
     applyModif(op.getPair(), op.getRouteIndex(), op.getPickupDeletion(), op.getDeliveryDeletion(), false);
 }
 
-
 bool CapacityConstraint::check(RemoveRoute const &op) const
 {
     return true;
-}   
+}
+
 void CapacityConstraint::apply(RemoveRoute const &op)
 {
     // add new empty route
@@ -148,16 +148,14 @@ void CapacityConstraint::print() const
 {
     std::cout << "Used capacity : \n";
     int i = 0;
-    for (const auto& capaVector : routeCapacities)
+    for (auto const &capaVector: routeCapacities)
     {
         std::cout << "#" << i << " : ";
-        for (const int capa : capaVector)
+        for (int const capa: capaVector)
         {
             std::cout << capa << ", ";
         }
-        std::cout << "\n ";
+        std::cout << "\n";
         i++;
     }
-    std::cout << "\n";
 }
-

src/lns/constraints/capacity/capacity_constraint.h

 #pragma once
-#include <vector>
-
 #include "lns/constraints/constraint.h"
-
-#include "lns/solution/solution.h"
 #include "lns/modification/pair/insert_pair.h"
 #include "lns/modification/pair/remove_pair.h"
 #include "lns/modification/route/insert_route.h"
 #include "lns/modification/route/remove_route.h"
+#include "lns/solution/solution.h"
+
+#include <vector>
 
 /**
  * Capacity constraint.
- * To check this, we keep track of the sum of all requests and check it against the max capacity.
- * TO DO, verification in O(1) by storing a matrix !
  */
 class CapacityConstraint : public Constraint
 {
 public:
     using CapacityVector = std::vector<int>;
-private:
 
+private:
+    /**
+     * Store the used capacity value when leaving one location in each route.
+     */
     std::vector<CapacityVector> routeCapacities;
 
     void apply(InsertPair const &op) override;
@@ -31,22 +31,21 @@ private:
     bool check(InsertRoute const &op) const override;
     bool check(RemovePair const &op) const override;
     bool check(RemoveRoute const &op) const override;
-public:
 
+public:
     explicit CapacityConstraint(Solution const &);
     CapacityConstraint(CapacityConstraint const &) = default;
 
     std::unique_ptr<Constraint> clone(Solution const &newOwningSolution) const override;
 
     /*
-    *   Given the solution, calculate the capacities when leaving each location
-    *   Modify the routeCapacities vector !
-    */
-    void initCapacities();
+     *  Compute from scratch the routeCapacities vector.
+     */
+    void initRouteCapacities();
 
     /*
-    *   Return capacity vector of a route
-    */
+     *  Return capacity vector of a route
+     */
     CapacityVector const &getRouteCapacities(int routeIndex) const;
 
     /*
@@ -62,7 +61,7 @@ public:
     /*
     *   Update the weight
     */
-    void applyModif(Pair const &pair, int routeIndex, int PickupPosition,  int DeliveryPosition, bool addPair);
+    void applyModif(Pair const &pair, int routeIndex, int PickupPosition, int DeliveryPosition, bool addPair);
 
     void print() const override;
 };
\ No newline at end of file

src/lns/constraints/time_window/forward_time_slack.cpp

+#include "forward_time_slack.h"
+
+#include "input/data.h"
+#include "lns/solution/solution.h"
+
+std::vector<TimeInteger> const &ForwardTimeSlack::getFTS() const
+{
+    return FTS;
+}
+
+std::vector<TimeInteger> const &ForwardTimeSlack::getEarliestArrival() const
+{
+    return earliestArrival;
+}
+
+std::vector<TimeInteger> const &ForwardTimeSlack::getLatestArrival() const
+{
+    return latestArrival;
+}
+
+void ForwardTimeSlack::initFTS(PDPTWData const &data, Route const &route)
+{
+    int n = route.getSize();
+    std::vector<int> const &routeIDs = route.getRoute();
+
+    earliestArrival = std::vector<TimeInteger>(n + 1, 0);
+    latestArrival = std::vector<TimeInteger>(n + 1, 0);
+    FTS = std::vector<TimeInteger>(n + 1, 0);
+
+    // Compute arrival times
+    double depotStart = data.getDepot().getTimeWindow().getStart();
+    double firstLocationStart = data.getLocation(routeIDs.at(0)).getTimeWindow().getStart();
+    double depotToFirst = data::TravelTime(data, 0, routeIDs.at(0));
+
+    earliestArrival.at(0) = std::max(firstLocationStart, depotStart + depotToFirst);
+
+    for (int i = 1; i < n; ++i)
+    {
+        double start = data.getLocation(routeIDs.at(i)).getTimeWindow().getStart();
+        double service = data.getLocation(routeIDs.at(i - 1)).getServiceDuration();
+        double travelTime = data::TravelTime(data, routeIDs.at(i - 1), routeIDs.at(i));
+
+        earliestArrival.at(i) = std::max(start, earliestArrival.at(i - 1) + service + travelTime);
+    }
+
+    // Compute earliest arrival at depot after last location
+    double lastService = data.getLocation(routeIDs.at(n - 1)).getServiceDuration();
+    double lastToDepot = data::TravelTime(data, routeIDs.at(n - 1), 0);
+    earliestArrival.at(n) = earliestArrival.at(n - 1) + lastService + lastToDepot;
+
+    // Compute latest arrival times
+    latestArrival.at(n - 1) = data.getLocation(routeIDs.at(n - 1)).getTimeWindow().getEnd();
+
+    for (int i = n - 2; i >= 0; --i)
+    {
+        double end = data.getLocation(routeIDs.at(i)).getTimeWindow().getEnd();
+        double service = data.getLocation(routeIDs.at(i)).getServiceDuration();
+        double travelTime = data::TravelTime(data, routeIDs.at(i), routeIDs.at(i + 1));
+
+        latestArrival.at(i) = std::min(latestArrival.at(i + 1) - service - travelTime, end);
+    }
+
+    // latest arrival at depot
+    latestArrival.at(n) = data.getDepot().getTimeWindow().getEnd();
+
+    // Compute FTS
+    for (int i = 0; i <= n; ++i)
+    {
+        FTS.at(i) = latestArrival.at(i) - earliestArrival.at(i);
+    }
+}
+
+bool ForwardTimeSlack::isPickupDeliveryInsertionValid(PDPTWData const &data, Route const &route, int pickupID,
+                                                      int deliveryID, int insertPickupIndex,
+                                                      int insertDeliveryIndex) const
+{
+    std::vector<int> const &routeIDs = route.getRoute();
+    int n = route.getSize();
+
+    // safety
+    if (insertPickupIndex > insertDeliveryIndex)
+    {
+        return false;
+    }
+
+    // if the route is empty, only check de time window
+    if (routeIDs.empty())
+    {
+        double startTW = data.getLocation(pickupID).getTimeWindow().getStart();
+        double endTW = data.getLocation(deliveryID).getTimeWindow().getEnd();
+
+        return (insertPickupIndex == 0) && (insertDeliveryIndex == 0) &&
+               (!(startTW > data.getDepot().getTimeWindow().getEnd() ||
+                  endTW < data.getDepot().getTimeWindow().getStart()));
+    }
+
+
+    // Algo
+
+    // Pickup
+    // two specific cases, insertion in first or last
+    int prev = (insertPickupIndex == 0) ? 0 : routeIDs.at(insertPickupIndex - 1);
+    int next = (insertPickupIndex == n) ? 0 : routeIDs.at(insertPickupIndex);
+
+    double arrivalPrev =
+            (prev == 0) ? data.getDepot().getTimeWindow().getStart() : getEarliestArrival().at(insertPickupIndex - 1);
+    double travelToNew = data::TravelTime(data, prev, pickupID);
+    double serviceNew = data.getLocation(pickupID).getServiceDuration();
+    double travelNewToNext = data::TravelTime(data, pickupID, next);
+
+    double newArrival = std::max(data.getLocation(pickupID).getTimeWindow().getStart(),
+                                 arrivalPrev + data.getLocation(prev).getServiceDuration() + travelToNew);
+
+
+    // Check TW
+    if (!(data.getLocation(pickupID).getTimeWindow().isValid(newArrival)))
+    {
+        return false;
+    }
+
+    // Calculate the pickup duration
+    double arrivalNext =
+            std::max(data.getLocation(next).getTimeWindow().getStart(), newArrival + serviceNew + travelNewToNext);
+
+    double pickupDuration = arrivalNext - earliestArrival.at(insertPickupIndex);
+    // Check FTS
+    if (/*next != 0 &&*/ pickupDuration >= getFTS().at(insertPickupIndex))
+    {
+        return false;
+    }
+
+    // Delivery
+    // if insertDeliveryIndex == insertPickupIndex, then Pickup is just before Delivery in the route
+    int prevDelivery = (insertDeliveryIndex == insertPickupIndex) ? pickupID : routeIDs.at(insertDeliveryIndex - 1);
+    int nextDelivery = (insertDeliveryIndex >= n) ? 0 : routeIDs.at(insertDeliveryIndex);
+
+    // We have to take into account the duration of the pickup
+    // (previous earliest arrival) + (time taken by the pickup) - (time of traveling between prev and next)
+    double arrivalPrevDelivery =
+            (prevDelivery == pickupID) ? newArrival : getEarliestArrival().at(insertDeliveryIndex - 1) + pickupDuration;
+
+    double travelToDelivery = data::TravelTime(data, prevDelivery, deliveryID);
+    double serviceDelivery = data.getLocation(deliveryID).getServiceDuration();
+    double travelDeliveryToNext = data::TravelTime(data, deliveryID, nextDelivery);
+    double newArrivalDelivery =
+            std::max(data.getLocation(deliveryID).getTimeWindow().getStart(),
+                     arrivalPrevDelivery + data.getLocation(prevDelivery).getServiceDuration() + travelToDelivery);
+
+    // Check TW
+    if (!(data.getLocation(deliveryID).getTimeWindow().isValid(newArrivalDelivery)))
+    {
+        return false;
+    }
+
+    // Calculate the delivery duration
+    double arrivalNextDelivery = std::max(data.getLocation(nextDelivery).getTimeWindow().getStart(),
+                                          newArrivalDelivery + serviceDelivery + travelDeliveryToNext);
+    double deliveryDuration = arrivalNextDelivery - earliestArrival.at(insertDeliveryIndex);
+
+    // Check FTS validity
+    if (/*nextDelivery != 0 &&*/ deliveryDuration >= getFTS().at(insertDeliveryIndex))
+    {
+        return false;
+    }
+
+    return true;
+}
+
+// TO DO some cuts in the loop to avoid recalculate everything
+void ForwardTimeSlack::updateFTSAfterInsertion(PDPTWData const &data, Route const &route, int insertPickupIndex,
+                                               int insertDeliveryIndex)
+{
+    int n = route.getSize();
+    std::vector<int> const &routeIDs = route.getRoute();
+
+    // Resize vectors to match the updated route size
+    earliestArrival.resize(n + 1, 0);
+    latestArrival.resize(n + 1, 0);
+    FTS.resize(n + 1, 0);
+
+    // Recalculate earliest arrival times
+    for (int i = insertPickupIndex; i <= n; ++i)
+    {
+        if (i == 0)// specific case for the first element
+        {
+            double depotStart = data.getDepot().getTimeWindow().getStart();
+            double firstLocationStart = data.getLocation(routeIDs.at(0)).getTimeWindow().getStart();
+            double depotToFirst = data::TravelTime(data, 0, routeIDs.at(0));
+            earliestArrival.at(0) = std::max(firstLocationStart, depotStart + depotToFirst);
+        }
+        else
+        {
+            int prev = routeIDs.at(i - 1);
+            int curr = (i == n) ? 0 : routeIDs.at(i);// specific case for the last element
+            double start = data.getLocation(curr).getTimeWindow().getStart();
+            double service = data.getLocation(prev).getServiceDuration();
+            double travelTime = data::TravelTime(data, prev, curr);
+
+            earliestArrival.at(i) = std::max(start, earliestArrival.at(i - 1) + service + travelTime);
+        }
+    }
+
+    // try add element in front then recalculate from the delivery index
+    // Recalculate latest arrival times
+    latestArrival.at(n) = data.getDepot().getTimeWindow().getEnd();
+    latestArrival.at(n - 1) = data.getLocation(routeIDs.at(n - 1)).getTimeWindow().getEnd();
+    for (int i = n - 2; i >= 0; --i)
+    {
+        double end = data.getLocation(routeIDs.at(i)).getTimeWindow().getEnd();
+        double service = data.getLocation(routeIDs.at(i)).getServiceDuration();
+        double travelTime = data::TravelTime(data, routeIDs.at(i), routeIDs.at(i + 1));
+
+        latestArrival.at(i) = std::min(latestArrival.at(i + 1) - service - travelTime, end);
+    }
+
+    // Recalculate FTS
+    for (int i = 0; i <= n; ++i)
+    {
+        FTS.at(i) = latestArrival.at(i) - earliestArrival.at(i);
+    }
+}
+
+// TO DO some cuts in the loop to avoid recalculate everything
+void ForwardTimeSlack::updateFTSAfterDeletion(PDPTWData const &data, Route const &route, int removePickupIndex,
+                                              int removeDeliveryIndex)
+{
+    int n = route.getSize();
+    std::vector<int> const &routeIDs = route.getRoute();
+
+    // empty route
+    if (n == 0)
+    {
+        earliestArrival.clear();
+        latestArrival.clear();
+        FTS.clear();
+        return;
+    }
+
+    // Resize vectors to match the updated route size
+    earliestArrival.resize(n + 1, 0);
+    latestArrival.resize(n + 1, 0);
+    FTS.resize(n + 1, 0);
+
+    // Recalculate earliest arrival times
+    for (int i = 0; i <= n; ++i)
+    {
+        if (i == 0)// Specific case for the first element
+        {
+            double depotStart = data.getDepot().getTimeWindow().getStart();
+            double firstLocationStart = data.getLocation(routeIDs.at(0)).getTimeWindow().getStart();
+            double depotToFirst = data::TravelTime(data, 0, routeIDs.at(0));
+            earliestArrival.at(0) = std::max(firstLocationStart, depotStart + depotToFirst);
+        }
+        else
+        {
+            int prev = routeIDs.at(i - 1);
+            int curr = (i == n) ? 0 : routeIDs.at(i);
+            double start = data.getLocation(curr).getTimeWindow().getStart();
+            double service = data.getLocation(prev).getServiceDuration();
+            double travelTime = data::TravelTime(data, prev, curr);
+
+            earliestArrival.at(i) = std::max(start, earliestArrival.at(i - 1) + service + travelTime);
+        }
+    }
+
+    // Recalculate latest arrival times
+    latestArrival.at(n) = data.getDepot().getTimeWindow().getEnd();
+    latestArrival.at(n - 1) = data.getLocation(routeIDs.at(n - 1)).getTimeWindow().getEnd();
+    for (int i = n - 2; i >= 0; --i)
+    {
+        double end = data.getLocation(routeIDs.at(i)).getTimeWindow().getEnd();
+        double service = data.getLocation(routeIDs.at(i)).getServiceDuration();
+        double travelTime = data::TravelTime(data, routeIDs.at(i), routeIDs.at(i + 1));
+
+        latestArrival.at(i) = std::min(latestArrival.at(i + 1) - service - travelTime, end);
+    }
+
+    // Recalculate FTS
+    for (int i = 0; i <= n; ++i)
+    {
+        FTS.at(i) = latestArrival.at(i) - earliestArrival.at(i);
+    }
+}
+
+void ForwardTimeSlack::print() const
+{
+    for (TimeInteger timeVal: earliestArrival)
+    {
+        std::cout << timeVal << " ";
+    }
+    std::cout << std::endl;
+    for (TimeInteger timeVal: latestArrival)
+    {
+        std::cout << timeVal << " ";
+    }
+    std::cout << std::endl;
+    for (TimeInteger ftsVal: FTS)
+    {
+        std::cout << ftsVal << " ";
+    }
+}
\ No newline at end of file

src/lns/constraints/time_window/forward_time_slack.h

+#pragma once
+
+#include "lns/solution/route.h"
+#include "types.h"
+
+/*
+ *  Represent the forward time slack of a route.
+ */
+class ForwardTimeSlack
+{
+private:
+    std::vector<TimeInteger> earliestArrival;
+    std::vector<TimeInteger> latestArrival;
+    std::vector<TimeInteger> FTS;
+
+public:
+    explicit ForwardTimeSlack() = default;
+
+    /*
+     *  Check if the insertion of a new pickup/delivery pair is feasible.
+     *  ForwardTimeSlack is supposed to be correct. Only check FTS and Time windows.
+     */
+    bool isPickupDeliveryInsertionValid(PDPTWData const &data, Route const &route, int pickupID, int deliveryID,
+                                        int insertPickupIndex, int insertDeliveryIndex) const;
+
+    /*
+     *  Compute FTS values from scratch.
+     */
+    void initFTS(PDPTWData const &data, Route const &route);
+
+    /*
+     *  Update FTS.
+     *  Suppose the route has already been updated (but not the FTS)
+     */
+    void updateFTSAfterInsertion(PDPTWData const &data, Route const &route, int insertPickupIndex,
+                                 int insertDeliveryIndex);
+    /*
+     *  Update FTS.
+     *  Suppose the route has already been updated (but not the FTS)
+     */
+    void updateFTSAfterDeletion(PDPTWData const &data, Route const &route, int removePickupIndex,
+        int removeDeliveryIndex);
+
+
+    std::vector<TimeInteger> const &getFTS() const;
+    std::vector<TimeInteger> const &getEarliestArrival() const;
+    std::vector<TimeInteger> const &getLatestArrival() const;
+
+    void print() const;
+};
\ No newline at end of file

src/lns/constraints/time_window/time_window_constraint.cpp

@@ -3,6 +3,7 @@
 #include "input/data.h"
 #include "input/pdptw_data.h"
 #include "input/time_window.h"
+#include "lns/constraints/time_window/forward_time_slack.h"
 #include "lns/modification/pair/insert_pair.h"
 #include "lns/modification/pair/remove_pair.h"
 #include "lns/modification/route/remove_route.h"
@@ -10,176 +11,71 @@
 
 TimeWindowConstraint::TimeWindowConstraint(Solution const &solution) : Constraint(solution)
 {
-    allRouteReachTimes.resize(getSolution().getRoutes().size());
+    FTSContainer.clear();
+    for (Route const &route: solution.getRoutes())
+    {
+        FTSContainer.emplace_back();
+    }
 }
 
 TimeWindowConstraint::~TimeWindowConstraint()
 {
-    allRouteReachTimes.clear();
+    //allRouteReachTimes.clear();
+    FTSContainer.clear();
 }
 
 std::unique_ptr<Constraint> TimeWindowConstraint::clone(Solution const &newOwningSolution) const
 {
     std::unique_ptr<TimeWindowConstraint> clonePtr = std::make_unique<TimeWindowConstraint>(newOwningSolution);
-    clonePtr->allRouteReachTimes = allRouteReachTimes;
+    //clonePtr->allRouteReachTimes = allRouteReachTimes;
+    clonePtr->FTSContainer = FTSContainer;
     return clonePtr;
 }
 
-void TimeWindowConstraint::computeReachTimes(PDPTWData const &data, std::vector<int> const &routeIDs,
-                                             ReachTimeVector &reachTimes) const
+void TimeWindowConstraint::initFTS()
 {
-    // Adjust the size of reachTimes vector
-    reachTimes.resize(routeIDs.size(), 0);
-    // Time to the first location
-    reachTimes.at(0) = data.getDepot().getTimeWindow().getStart() + data::TravelTime(data, 0, routeIDs.at(0));
-    // Compute other reachTimes (max between arrival and start of the time window + previous service time + travel time)
-    for (int i = 1; i < routeIDs.size(); ++i)
-    {
-        TimeInteger travelTime = data::TravelTime(data, routeIDs.at(i - 1), routeIDs.at(i));
-        // locations are indexed from 0 to n-1,
-        TimeInteger serviceTime = data.getLocation(routeIDs.at(i - 1)).getServiceDuration();
-        TimeInteger startTW = data.getLocation(routeIDs.at(i - 1)).getTimeWindow().getStart();
-
-        reachTimes.at(i) = std::max(reachTimes.at(i - 1), startTW) + serviceTime + travelTime;
-    }
-}
-
-void TimeWindowConstraint::initReachTimes()
-{
-    allRouteReachTimes = std::vector<ReachTimeVector>();
+    FTSContainer = std::vector<ForwardTimeSlack>();
+    //reserve
     int i = 0;
     for (Route const &route: getSolution().getRoutes())
     {
-        // init and calculate reach time value
-        allRouteReachTimes.emplace_back();
-        computeReachTimes(getSolution().getData(), route.getRoute(), allRouteReachTimes.at(i));
+        FTSContainer.emplace_back();
+        FTSContainer.at(i).initFTS(getSolution().getData(), route);
         ++i;
     }
 }
 
-// Copie le vecteur de reach time de la route concernée
-// Insère les positions pickup/delivery
-// refait l'ordo sur le nouveau vecteur
 bool TimeWindowConstraint::checkInsertion(PDPTWData const &data, Pair const &pair, int routeIndex, int pickupPos,
                                           int deliveryPos) const
 {
-    ReachTimeVector const &reachTimes = allRouteReachTimes.at(routeIndex);
-
-    // COPY route vector
-    std::vector<int> route(getSolution().getRoute(routeIndex).getRoute().begin(),
-                           getSolution().getRoute(routeIndex).getRoute().end());
-    // COPY reachTimes vector
-    ReachTimeVector newReachTimes(reachTimes.begin(), reachTimes.end());
-
-    // Insert pickup and delivery
-    route.insert(route.begin() + deliveryPos, pair.getDelivery().getId());
-    route.insert(route.begin() + pickupPos, pair.getPickup().getId());
-
-
-    // std::cout << "\n";
-    // for (auto pos : route)
-    // {
-    //     std::cout << pos << " ";
-    // }
-    // std::cout << "\n";
-
-    // Compute new reach time
-    computeReachTimes(data, route, newReachTimes);
-
-    // std::cout << "\n";
-    // for (auto pos : newReachTimes)
-    // {
-    //     std::cout << pos << " ";
-    // }
-    // std::cout << "\n";
-
-    // Check Time Windows
-    for (int i = 0; i < newReachTimes.size(); ++i)
-    {
-        if (!data.getLocation(route.at(i)).getTimeWindow().isValid(newReachTimes.at(i)))
-        {
-            return false;
-        }
-    }
-    return true;
+    return FTSContainer.at(routeIndex)
+            .isPickupDeliveryInsertionValid(data,
+                                            getSolution().getRoute(routeIndex),
+                                            pair.getPickup().getId(),
+                                            pair.getDelivery().getId(),
+                                            pickupPos,
+                                            deliveryPos);
 }
 
-// A lot of check here because i had problems with this function
+
 void TimeWindowConstraint::ApplyModif(PDPTWData const &data, Pair const &pair, int routeIndex, int pickupPos,
                                       int deliveryPos, bool addPair)
 {
-
-    // Check the routeIndex validity
-    if (routeIndex < 0 || routeIndex >= getSolution().getRoutes().size())
-    {
-        spdlog::error("Error: routeIndex out of bounds ({})", routeIndex);
-        return;
-    }
-
-    // Copy of the route vector (problem that cause a duplication of the TW)
-    std::vector<int> routeIDs = getSolution().getRoute(routeIndex).getRoute();
-
-    // // Check Position validity
-    // if (pickupPos < 0 || pickupPos > routeIDs.size() ||
-    //     deliveryPos < 0 || deliveryPos > routeIDs.size()) {
-    //     spdlog::error("Error: Indices pickupPos ({}) or deliveryPos ({}) are invalid for a route size of {}.",
-    //         pickupPos, deliveryPos, routeIDs.size());
-    //     return;
-    // }
-
-    // if (addPair) {
-    //     routeIDs.insert(routeIDs.begin() + deliveryPos, pair.getDelivery().getId());
-    //     routeIDs.insert(routeIDs.begin() + pickupPos, pair.getPickup().getId());
-    // }
-    // else {
-    //     if (deliveryPos < routeIDs.size() && pickupPos < routeIDs.size() && pickupPos != deliveryPos) {
-    //         if (deliveryPos > pickupPos) {
-    //             routeIDs.erase(routeIDs.begin() + deliveryPos);
-    //             routeIDs.erase(routeIDs.begin() + pickupPos);
-    //         } else {
-    //             routeIDs.erase(routeIDs.begin() + pickupPos);
-    //             routeIDs.erase(routeIDs.begin() + deliveryPos);
-    //         }
-    //     } else {
-    //         spdlog::error("Error: Invalid indices for removal (pickupPos: {}, deliveryPos: {}).", pickupPos, deliveryPos);
-    //         return;
-    //     }
-    // }
-
-    // the route is empty !
-    if (routeIDs.empty())
+    if (addPair)
     {
-        allRouteReachTimes.at(routeIndex).clear();
+        FTSContainer.at(routeIndex)
+                .updateFTSAfterInsertion(data, getSolution().getRoute(routeIndex), pickupPos, deliveryPos);
     }
     else
     {
-        allRouteReachTimes.at(routeIndex).resize(routeIDs.size(), 0);
-
-        allRouteReachTimes.at(routeIndex).at(0) =
-                data.getDepot().getTimeWindow().getStart() + data::TravelTime(data, 0, routeIDs.at(0));
-
-        // Compute reach times
-        for (size_t i = 1; i < routeIDs.size(); ++i)
-        {
-            if (i - 1 >= routeIDs.size() || i >= routeIDs.size())
-            {
-                spdlog::error("Error: Out-of-bounds access to routeIDs in the computation loop.");
-                return;
-            }
-
-            TimeInteger travelTime = data::TravelTime(data, routeIDs.at(i - 1), routeIDs.at(i));
-            TimeInteger serviceTime = data.getLocation(routeIDs.at(i - 1)).getServiceDuration();
-            TimeInteger startTW = data.getLocation(routeIDs.at(i - 1)).getTimeWindow().getStart();
-
-            allRouteReachTimes.at(routeIndex).at(i) =
-                    std::max(allRouteReachTimes.at(routeIndex).at(i - 1), startTW) + serviceTime + travelTime;
-        }
+        FTSContainer.at(routeIndex)
+                .updateFTSAfterDeletion(data, getSolution().getRoute(routeIndex), pickupPos, deliveryPos);
     }
+
 }
 
 bool TimeWindowConstraint::check(InsertPair const &op) const
 {
-    //std::cout << " #TW Check";
     return checkInsertion(getSolution().getData(),
                           op.getPair(),
                           op.getRouteIndex(),
@@ -189,7 +85,6 @@ bool TimeWindowConstraint::check(InsertPair const &op) const
 
 void TimeWindowConstraint::apply(InsertPair const &op)
 {
-    //std::cout << "-> Apply Modification on Time Windows \n";
     ApplyModif(getSolution().getData(),
                op.getPair(),
                op.getRouteIndex(),
@@ -205,7 +100,7 @@ bool TimeWindowConstraint::check(InsertRoute const &op) const
 
 void TimeWindowConstraint::apply(InsertRoute const &op)
 {
-    allRouteReachTimes.emplace_back();
+    FTSContainer.emplace_back();
 }
 
 bool TimeWindowConstraint::check(RemovePair const &op) const
@@ -230,34 +125,18 @@ bool TimeWindowConstraint::check(RemoveRoute const &op) const
 
 void TimeWindowConstraint::apply(RemoveRoute const &op)
 {
-    allRouteReachTimes.erase(allRouteReachTimes.begin() + op.getRouteIndex());
-}
-
-std::vector<TimeWindowConstraint::ReachTimeVector> const &TimeWindowConstraint::getallRouteReachTimes() const
-{
-    return allRouteReachTimes;
+    FTSContainer.erase(FTSContainer.begin() + op.getRouteIndex());
 }
 
 void TimeWindowConstraint::print() const
 {
-    std::cout << "Reach Times : \n";
+    std::cout << "Time Window : Earliest / Latest / FTS" << std::endl;
     int i = 0;
-    for (auto const &reachTimes: getallRouteReachTimes())
+    for (ForwardTimeSlack FTS: FTSContainer)
     {
-        //std::cout << reachTimes.size() << ", ";
-        std::cout << "#" << i << " : ";
-        for (const TimeInteger reachTime: reachTimes)
-        {
-            std::cout << reachTime << ", ";
-        }
-        std::cout << "\n";
-        std::cout << "#" << i << " : ";
-        for (int locID: getSolution().getRoute(i).getRoute())
-        {
-            std::cout << locID << ", ";
-        }
-        std::cout << "\n ";
+        std::cout << "#" << i << std::endl;
+        FTS.print();
+        std::cout << std::endl;
         i++;
     }
-    std::cout << "\n";
 }

src/lns/constraints/time_window/time_window_constraint.h

 #pragma once
 
-#include "lns/constraints/constraint.h"
 #include "input/pair.h"
 #include "input/time_window.h"
+#include "lns/constraints/constraint.h"
+#include "lns/constraints/time_window/forward_time_slack.h"
 #include "lns/solution/route.h"
+
 #include <vector>
 
 /**
@@ -17,59 +19,44 @@ class TimeWindowConstraint : public Constraint
 public:
     using ReachTimeVector = std::vector<TimeInteger>;
 
+
     explicit TimeWindowConstraint(Solution const &);
     TimeWindowConstraint(TimeWindowConstraint const &) = default;
     ~TimeWindowConstraint() override;
 
-    // Given the solution, calculate the reachTimes of each location.
-    // Modify the routeCapacities vector !
-    void initReachTimes();
+
+    /**
+     *  Initialize the FTS vectors
+     */
+    void initFTS();
 
     std::unique_ptr<Constraint> clone(Solution const &newOwningSolution) const override;
-    
-    const std::vector<ReachTimeVector> & getallRouteReachTimes() const;
 
-    void print() const override;
-private:
 
-    /**
-    *   For each route, store the reach time of each location
-    */
-    std::vector<ReachTimeVector> allRouteReachTimes;  
-
-    // Recalcul tous les reachTimes pour une route
-    ReachTimeVector computeReachTimes(const Route& route) const;
-    bool updateReachTimes(const PDPTWData& data, Pair & pair,int routeIndex, int pickupPos, int deliveryPos);
+    void print() const override;
 
-    /**
-    *   Check if the insertion of a pair pickup/delivery is valid or not.
-    *   COPY the route where we insert the pair !
-    */
-    bool checkInsertion(const PDPTWData& data, const Pair & pair, int routeIndex, int pickupPos, int deliveryPos) const; 
+private:
+    std::vector<ForwardTimeSlack> FTSContainer;
 
     /**
-    *   Given route IDs, compute the reachTimes of the route.
-    *   If a vehicule is in advance, set the reach time to start of the timeWindow.
-    *   Modify reachTimes vector !
-    */
-    void computeReachTimes(const PDPTWData& data, const std::vector<int> & routeIDs, ReachTimeVector & reachTimes) const;
+     *  Check if the insertion of a pair pickup/delivery is valid or not.
+     *  COPY the route where we insert the pair !
+     */
+    bool checkInsertion(PDPTWData const &data, Pair const &pair, int routeIndex, int pickupPos, int deliveryPos) const;
 
     /**
-    *   Given route IDs, modify the reachTimes
-    */
-    void ApplyModif(const PDPTWData& data, const Pair & pair, int routeIndex, int pickupPos, int deliveryPos,  bool addPair);
-
+     *  Apply the modification to the time window constraint
+     */
+    void ApplyModif(PDPTWData const &data, Pair const &pair, int routeIndex, int pickupPos, int deliveryPos,
+                    bool addPair);
 
     void apply(InsertPair const &op) override;
     void apply(InsertRoute const &op) override;
     void apply(RemovePair const &op) override;
     void apply(RemoveRoute const &op) override;
- 
+
     bool check(InsertPair const &op) const override;
     bool check(InsertRoute const &op) const override;
     bool check(RemovePair const &op) const override;
     bool check(RemoveRoute const &op) const override;
-
-public: 
-
 };
\ No newline at end of file

src/lns/lns.cpp

@@ -4,6 +4,8 @@
 #include "lns/operators/selector/operator_selector.h"
 #include "output/solution_checker.h"
 
+#include <chrono>
+
 namespace
 {
     bool isBetterSolution(Solution const &candidateSolution, Solution const &bestKnownSol)
@@ -11,22 +13,82 @@ namespace
         return candidateSolution.getCost() < bestKnownSol.getCost();
     }
 
+    using lns_clock = std::chrono::high_resolution_clock;
+    using lns_time_point = std::chrono::time_point<lns_clock, std::chrono::nanoseconds>;
+
+    struct LnsRuntimeData
+    {
+        Solution bestSolution;
+        unsigned int numberOfIteration = 0;
+        lns_time_point start = lns_clock::now();
+    };
+
+    /**
+     * @return the number of seconds between point and now (the moment the function is called).
+     */
+    unsigned long getTimeSinceInSec(lns_time_point point)
+    {
+        return std::chrono::duration_cast<std::chrono::seconds>(lns_clock::now() - point).count();
+    }
+
+    unsigned long getTimeSinceInMs(lns_time_point point)
+    {
+        return std::chrono::duration_cast<std::chrono::milliseconds>(lns_clock::now() - point).count();
+    }
+
+    bool triggerLogProgress(LnsRuntimeData const &)
+    {
+        static auto lastTrigger = lns_clock::now();
+        if (getTimeSinceInSec(lastTrigger) >= 10)
+        {
+            lastTrigger = lns_clock::now();
+            return true;
+        }
+        return false;
+    }
+
+    void logProgress(LnsRuntimeData const &runtime, Solution const &actualSolution)
+    {
+        if (triggerLogProgress(runtime))
+        {
+            unsigned long actualTime = getTimeSinceInMs(runtime.start);
+            double iterPerSecond = 1000 * runtime.numberOfIteration / static_cast<double>(actualTime);
+            std::string speedLog;
+            if (iterPerSecond < 1)
+            {
+                speedLog = fmt::format("{}s/100 iterations", 100 / iterPerSecond);
+            }
+            else
+            {
+                speedLog = fmt::format("{:.1f} i/s", iterPerSecond);
+            }
+            std::string missingRequestLog;
+            long requestsMissing = runtime.bestSolution.missingPairCount();
+            //... log TO DO
+        }
+    }
+
 }// namespace
 
 output::LnsOutput lns::runLns(Solution const &initialSolution, OperatorSelector &opSelector,
-                         AcceptanceFunction const &acceptFunctor)
+                              AcceptanceFunction const &acceptFunctor)
 {
     /**
      * The solution is the base solution used to create the neighbor solution.
      * It is constant unless we accept the candidate solution.
      */
     Solution actualSolution = initialSolution;
-    Solution bestSolution = initialSolution;
+    LnsRuntimeData runtime = {actualSolution};
+
+    // temporary fixed iteration
+    int iterationMax = 100;
+    while (iterationMax > 0)
+    {
+        // Init iteration
+        ++runtime.numberOfIteration;
+        logProgress(runtime, actualSolution);
 
-    int it = 100;
-    while (it > 0)
-    {   
-        std::cout << "\n" << 101-it << " : ";
+        std::cout << "\n" << runtime.numberOfIteration << " : ";
         /**
          * The solution on which we apply the operators.
          * It is discarded at the end of each loop if it is not accepted by the Acceptance Function.
@@ -38,26 +100,28 @@ output::LnsOutput lns::runLns(Solution const &initialSolution, OperatorSelector
         // Apply operators
         destructReconstructPair.destructor().destroySolution(candidateSolution);
         destructReconstructPair.reconstructor().reconstructSolution(candidateSolution, 0.01);
-
+        candidateSolution.computeAndStoreSolutionCost();
         // Update best solution
-        if (isBetterSolution(candidateSolution, bestSolution))
+        if (isBetterSolution(candidateSolution, runtime.bestSolution))
         {
             std::cout << "\n > new Best Solution \n";
             checker::checkAll(candidateSolution, candidateSolution.getData(), false);
 
-            bestSolution = candidateSolution;
+            runtime.bestSolution = candidateSolution;
             opSelector.betterSolutionFound();
         }
 
         // Check if we use the candidate solution as the new actual solution
         // operator can force to take the new solution
         if (destructReconstructPair.forceTakeSolution() ||
-            acceptFunctor(candidateSolution, actualSolution, bestSolution) == AcceptationStatus::ACCEPT)
+            acceptFunctor(candidateSolution, actualSolution, runtime.bestSolution) == AcceptationStatus::ACCEPT)
         {
             actualSolution = std::move(candidateSolution);
         }
-        --it;
+        --iterationMax;
     }
-    auto result = output::LnsOutput(std::move(bestSolution), it, it);
+
+    auto result = output::LnsOutput(
+            std::move(runtime.bestSolution), runtime.numberOfIteration, getTimeSinceInSec(runtime.start));
     return result;
 }
\ No newline at end of file

src/lns/modification/pair/remove_pair.cpp

 #include "remove_pair.h"
 
 #include "input/data.h"
+#include "lns/solution/solution.h"
 
 RemovePair::RemovePair(int routeIndex, int pickupDeletion, int deliveryDeletion, Pair const &pair)
     : routeIndex(routeIndex), pickupDeletion(pickupDeletion), deliveryDeletion(deliveryDeletion),

src/lns/modification/route/remove_route.cpp

@@ -18,16 +18,19 @@ RemoveRoute::RemoveRoute(int routeIndex, std::vector<int> &&removedPairID)
 void RemoveRoute::modifySolution(Solution &solution)
 {
     std::vector<Route> &routes = solution.getRoutes();
-    std::vector<int> const &locationIDs = routes.at(routeIndex).getRoute();
+    if (routes.at(routeIndex).getRoute().empty())
+    {
+        std::vector<int> const &locationIDs = routes.at(routeIndex).getRoute();
 
-    // update removedPairID
-    removedPairID.reserve(routes.at(routeIndex).getSize());
+        // update removedPairID
+        removedPairID.reserve(routes.at(routeIndex).getSize());
 
-    for (int id : locationIDs)
-    {
-        if (solution.getData().getLocation(id).getLocType() == LocType::PICKUP)
+        for (int id: locationIDs)
         {
-            removedPairID.push_back(id);
+            if (solution.getData().getLocation(id).getLocType() == LocType::PICKUP)
+            {
+                removedPairID.push_back(id);
+            }
         }
     }
     routes.erase(routes.begin() + routeIndex);
@@ -46,4 +49,9 @@ int RemoveRoute::getRouteIndex() const
 std::vector<int> const &RemoveRoute::getDeletedPairs() const
 {
     return removedPairID;
+}
+
+ModificationApplyVariant RemoveRoute::asApplyVariant() const
+{
+    return *this;
 }
\ No newline at end of file

src/lns/modification/route/remove_route.h

@@ -36,5 +36,8 @@ public:
     */
     std::vector<int> const &getDeletedPairs() const override;
 
+    ModificationApplyVariant asApplyVariant() const override;
+    
+
     int getRouteIndex() const;
 };
\ No newline at end of file

src/lns/operators/destruction/clean_empty_route.cpp

+#include "clean_empty_route.h"
+#include "lns/modification/route/remove_route.h"
+
+CleanEmptyRoute::CleanEmptyRoute() {}
+
+void CleanEmptyRoute::destroySolution(Solution &solution) const
+{
+    std::cout << "Clean ";
+    for (int routeIndex = solution.getRoutes().size() - 1; routeIndex >= 0; --routeIndex)
+    {
+        if (solution.getRoute(routeIndex).getRoute().empty())
+        {
+            RemoveRoute remRoute = RemoveRoute(routeIndex);
+            solution.applyDestructSolution(remRoute);
+        }
+    }
+}
\ No newline at end of file

src/lns/operators/destruction/clean_empty_route.h

+#pragma once
+
+#include "lns/operators/abstract_operator.h"
+
+class CleanEmptyRoute : public DestructionOperator
+{
+public:
+    explicit CleanEmptyRoute();
+
+    /**
+     * This operator removes empty routes from the solution.
+     */
+    void destroySolution(Solution &solution) const override;
+
+};

src/lns/operators/destruction/random_destroy.cpp

@@ -18,7 +18,7 @@ void RandomDestroy::destroySolution(Solution &solution) const
     while (remainingPairToDelete > 0)
     {
         // too complicated
-        // Other (simpler) option -> choose random route then choose random pair ?
+        // Other (simpler) option -> choose random route then choose random pair ? (does not respect equiprobability)
      
         // choose a random location
         int locationNumber = util::getRandomInt(0, remainingPairToDelete * 2 - 1);

src/lns/operators/destruction/string_removal.cpp

+#include "string_removal.h"
+
+#include "input/location.h"
+#include "input/pair.h"
+#include "input/pdptw_data.h"
+#include "lns/modification/pair/remove_pair.h"
+#include "utils.h"
+
+#include <vector>
+
+namespace
+{
+
+    double computeAverageCardinality(std::vector<Route> const &routes)
+    {
+        if (routes.empty())
+        {
+            return 0.0;
+        }
+        return std::accumulate(routes.begin(),
+                               routes.end(),
+                               0.0,
+                               [](double sum, Route const &route) { return sum + route.getSize(); }) /
+               routes.size();
+    }
+
+    /**
+     * Remove a string in the solution
+     * For now it uses startLocation as the center of the string
+     * @param routeIndex the path on which the removal happens
+     * @param stringLength the number of requests to remove, expect it to be inferior to the number of requests in the path
+     * @param startLocation the location found to start the string, does not mean it will the smallest index removed
+     */
+    void removeString(Solution &solution, int routeIndex, unsigned int stringLength, int startLocation)
+    {
+        // ID of pair to remove
+        std::vector<int> toRemove;
+        toRemove.reserve(stringLength);
+
+        PDPTWData const &data = solution.getData();
+        Route const &route = solution.getRoute(routeIndex);
+        std::vector<int> const &routeLocationIDs = route.getRoute();
+
+        // Position in the route of the first location
+        int startPosition = route.getIndex(startLocation);
+
+        // the location is not in a route
+        if (startPosition == -1)
+        {
+            return;
+        }
+
+        int left = std::max(0, startPosition - ((int) stringLength / 2));
+        int right = left + (int) stringLength - 1;
+
+        // we don't want to go out of bound
+        if (right >= routeLocationIDs.size())
+        {
+            right = routeLocationIDs.size() - 1;
+            left = std::max(0, right - (int) stringLength + 1);// keep a number of stringLength elements
+        }
+
+        // collect each locationID and store the pairID in the toRemove vector
+        for (int i = left; i <= right; ++i)
+        {
+            // pickupID = pairID
+            // if pickup and not already in the list
+            if (data.getLocation(routeLocationIDs.at(i)).getLocType() == LocType::PICKUP &&
+                std::find(toRemove.begin(), toRemove.end(), routeLocationIDs.at(i)) == toRemove.end())
+            {
+                toRemove.emplace_back(routeLocationIDs.at(i));
+            }
+            // or if delivery and pairID not already in the list
+            else if (data.getLocation(routeLocationIDs.at(i)).getLocType() == LocType::DELIVERY &&
+                     std::find(toRemove.begin(), toRemove.end(), data.getLocation(routeLocationIDs.at(i)).getPair()) ==
+                             toRemove.end())
+            {
+                toRemove.emplace_back(data.getLocation(routeLocationIDs.at(i)).getPair());
+            }
+        }
+
+        //remove the collected pairID
+        for (int pairID: toRemove)
+        {
+            // get index and Pair
+            int position = solution.getRoute(routeIndex).getIndex(pairID);
+            Index index = std::make_tuple(routeIndex, position, route.getPairLocationPosition(position, data));
+            Pair const &pair = data.getPair(pairID);
+
+            // apply remove
+            RemovePair removePair(index, pair);
+            solution.applyDestructSolution(removePair);
+        }
+    }
+
+    // Procedure RUIN(s, Lmax, c^bar)
+    //   1. Compute lmax_s, kmax_s, k_s according to equations (5), (6), and (7)
+    //   2. cseed_s ← RandomClient(s)
+    //   3. R ← ∅  // Set of affected tours
+    //   4. For each client c in Adjacent(cseed_s) Do
+    //       5. If c ∉ A and |R| < k_s Then
+    //           6. c*_t ← c
+    //           7. Compute lmax_t, lt according to equations (8) and (9)
+    //           8. A ← A ∪ RemoveClient(T, lt, c*_t)
+    //           9. R ← R ∪ {t}
+    //       10. End If
+    //   11. End For
+    // End Procedure
+
+    /**
+     *  Ruin procedure from the SISR paper.
+     *  c^bar = averageNumberRemovedElement.
+     *  L^max = maxStringSize.
+     */
+    void SISRsRuin(Solution &solution, unsigned int maxStringSize, unsigned int averageNumberRemovedElement)
+    {
+        // |t∈T| average number of location by route
+        auto averageRouteCardinality = static_cast<unsigned int>(computeAverageCardinality(solution.getRoutes()));
+
+        unsigned int maxSizeOfString = std::min(maxStringSize, averageRouteCardinality);               // (5) lmax_s
+        unsigned int maxNumberOfString = (4 * averageNumberRemovedElement) / (1 + maxSizeOfString) - 1;// (6) kmax_s
+        unsigned int numberOfString = util::getRandomInt(1, maxNumberOfString + 1);                    // (7) k_s
+
+        int locationSeed = util::getRandomInt(1, solution.getData().getLocationCount());
+
+        // store the index of routes where a string was removed
+        std::vector<int> routeIndexUsed;
+        routeIndexUsed.reserve(numberOfString);
+
+        // getClosestLocationsID returns the list of all location sorted from the closest to furthest
+        for (int neighbor: solution.getData().getClosestLocationsID(locationSeed))
+        {
+            // recover the routeIndex associated to the locationID, returns -1 if in the bank
+            int routeIndex = solution.getRouteIDOf(neighbor);
+
+            if (routeIndex != -1 && std::ranges::find(routeIndexUsed, routeIndex) == routeIndexUsed.end())
+            {
+                // (8) lmax_t
+                unsigned int maxSizeOfThisString =
+                        std::min(static_cast<int>(maxSizeOfString), solution.getRoute(routeIndex).getSize());
+                // (9) l_t
+                unsigned int actualSizeOfThisString = util::getRandomInt(1, maxSizeOfThisString);
+
+                removeString(solution, routeIndex, actualSizeOfThisString, neighbor);
+                routeIndexUsed.emplace_back(routeIndex);
+
+                if (routeIndexUsed.size() >= numberOfString)
+                {
+                    break;
+                }
+            }
+        }
+    }
+
+}// namespace
+
+void StringRemoval::destroySolution(Solution &solution) const
+{
+    std::cout << "SISR ";
+    SISRsRuin(solution, maxCardinalityOfString, averageNumberRemovedElement);
+}

src/lns/operators/destruction/string_removal.h

+#pragma once
+
+#include "lns/operators/abstract_operator.h"
+
+class StringRemoval : public DestructionOperator
+{
+public:
+    explicit StringRemoval(unsigned int maxCardinalityOfString, unsigned int averageNumberRemovedElement)
+        : maxCardinalityOfString(maxCardinalityOfString), averageNumberRemovedElement(averageNumberRemovedElement)
+    {}
+
+    // Each string has a max number of requests: MAX_CARDINALITY_OF_STRING
+    // L_max in the paper.
+    unsigned int maxCardinalityOfString = 10;
+    // This operator will remove in average AVERAGE_CUSTOMER_REMOVAL of requests
+    // c with a bar on top in the paper
+    unsigned int averageNumberRemovedElement = 10;
+
+    /**
+     * This operator removes numberOfPairsToDestroy pairs randomly in the solution.
+     */
+    void destroySolution(Solution &solution) const override;
+
+private:
+    int numberOfPairsToDestroy;
+};