From 1557e142edac9230246a3a917a82e4b34ff8d68b Mon Sep 17 00:00:00 2001
From: Laio Oriel Seman <laioseman@gmail.com>
Date: Fri, 27 Sep 2024 20:19:33 -0300
Subject: [PATCH] n-dimensional splay tree for fast bucket index computation

---
 bucket/BucketGraph.h              |  54 +++--
 bucket/BucketSolve.h              |  38 +++-
 bucket/BucketUtils.h              | 162 ++++++--------
 examples/VRPTW.cpp                |   4 +-
 examples/VRPTW.h                  |   2 +-
 external/hgs_vrptw/Individual.cpp |   4 +-
 external/hgs_vrptw/Individual.h   |  15 +-
 external/hgs_vrptw/Population.cpp |   2 +-
 include/Definitions.h             | 350 ++++++++++++++++++++++++++++++
 include/HGS.h                     |  18 +-
 src/BucketGraph.cpp               |  15 +-
 11 files changed, 519 insertions(+), 145 deletions(-)

diff --git a/bucket/BucketGraph.h b/bucket/BucketGraph.h
index 4b3c48f7..27b7ddb5 100644
--- a/bucket/BucketGraph.h
+++ b/bucket/BucketGraph.h
@@ -267,6 +267,8 @@ class BucketGraph {
 
     std::thread rih_thread;
     std::mutex  mtx; // For thread-safe access to merged_labels_improved
+    //
+    int redefine_counter = 0;
 
     std::vector<std::vector<int>> fw_ordered_sccs;
     std::vector<std::vector<int>> bw_ordered_sccs;
@@ -338,6 +340,13 @@ class BucketGraph {
     double min_red_cost = std::numeric_limits<double>::infinity();
     bool   first_reset  = true;
 
+    std::vector<int> dominance_checks_per_bucket;
+    int              non_dominated_labels_per_bucket;
+
+    // Interval tree to store bucket intervals
+    std::unordered_map<int, SplayTree> fw_job_interval_trees;
+    std::unordered_map<int, SplayTree> bw_job_interval_trees;
+
 #ifdef SCHRODINGER
     /**
      * @brief Retrieves a list of paths with negative reduced costs.
@@ -495,30 +504,30 @@ class BucketGraph {
      *
      * @param bucketInterval The new interval for the buckets.
      */
-    [[maybe_unused]] void redefine(int bucketInterval) {
+    void redefine(int bucketInterval) {
         this->bucket_interval = bucketInterval;
         intervals.clear();
         for (int i = 0; i < R_SIZE; ++i) { intervals.push_back(Interval(bucketInterval, 0)); }
 
-        define_buckets<Direction::Forward>();
-        define_buckets<Direction::Backward>();
-
-        fixed_arcs.resize(getJobs().size());
-        for (int i = 0; i < getJobs().size(); ++i) { fixed_arcs[i].resize(getJobs().size()); }
-
-        // make every fixed_buckets also have size buckets.size()
-        fw_fixed_buckets.resize(fw_buckets.size());
-        bw_fixed_buckets.resize(fw_buckets.size());
-
-        for (auto &fb : fw_fixed_buckets) { fb.resize(fw_buckets.size()); }
-        for (auto &bb : bw_fixed_buckets) { bb.resize(bw_buckets.size()); }
-        // set fixed_buckets to 0
-        for (auto &fb : fw_fixed_buckets) {
-            for (std::size_t i = 0; i < fb.size(); ++i) { fb[i] = 0; }
-        }
-        for (auto &bb : bw_fixed_buckets) {
-            for (std::size_t i = 0; i < bb.size(); ++i) { bb[i] = 0; }
-        }
+        reset_fixed();
+        reset_fixed_buckets();
+
+        PARALLEL_SECTIONS(
+            bi_sched,
+            SECTION {
+                // Section 1: Forward direction
+                define_buckets<Direction::Forward>();
+                fw_fixed_buckets.clear();
+                fw_fixed_buckets.resize(fw_buckets_size);
+                for (auto &fb : fw_fixed_buckets) { fb.assign(fw_buckets_size, 0); }
+            },
+            SECTION {
+                // Section 2: Backward direction
+                define_buckets<Direction::Backward>();
+                bw_fixed_buckets.clear();
+                bw_fixed_buckets.resize(bw_buckets_size);
+                for (auto &bb : bw_fixed_buckets) { bb.assign(bw_buckets_size, 0); }
+            });
 
         generate_arcs();
         for (auto &VRPJob : jobs) { VRPJob.sort_arcs(); }
@@ -575,6 +584,11 @@ class BucketGraph {
         for (auto &row : fixed_arcs) { std::fill(row.begin(), row.end(), 0); }
     }
 
+    void reset_fixed_buckets() {
+        for (auto &fb : fw_fixed_buckets) { std::fill(fb.begin(), fb.end(), 0); }
+        for (auto &bb : bw_fixed_buckets) { std::fill(bb.begin(), bb.end(), 0); }
+    }
+
     /**
      * @brief Checks the feasibility of a given forward and backward label.
      *
diff --git a/bucket/BucketSolve.h b/bucket/BucketSolve.h
index bea4d257..8b18efaa 100644
--- a/bucket/BucketSolve.h
+++ b/bucket/BucketSolve.h
@@ -219,6 +219,12 @@ std::vector<double> BucketGraph::labeling_algorithm() noexcept {
                             const auto &to_bucket_labels =
                                 buckets[to_bucket].get_labels(); // Get existing labels in the destination bucket
 
+                            if constexpr (S == Stage::Four) {
+                                // Track dominance checks for this bucket
+                                if constexpr (D == Direction::Forward) {
+                                    dominance_checks_per_bucket[to_bucket] += to_bucket_labels.size();
+                                }
+                            }
                             // Stage-specific dominance check
                             if constexpr (S == Stage::One) {
                                 // If the new label has lower cost, remove dominated labels
@@ -414,7 +420,10 @@ std::vector<Label *> BucketGraph::bi_labeling_algorithm() {
     for (auto bucket = 0; bucket < fw_buckets_size; ++bucket) {
         auto       &current_bucket = fw_buckets[bucket];          // Get the current bucket
         const auto &labels         = current_bucket.get_labels(); // Get labels in the current bucket
-
+        //
+        if constexpr (S == Stage::Four) {
+            non_dominated_labels_per_bucket += labels.size(); // Track non-dominated labels
+        }
         // Process each label in the bucket
         for (const Label *L : labels) {
             // if (L->resources[TIME_INDEX] > q_star[TIME_INDEX]) { continue; } // Skip if label exceeds q_star
@@ -483,7 +492,30 @@ std::vector<Label *> BucketGraph::bi_labeling_algorithm() {
         sPool.iterate();
     }
 #endif
+    if constexpr (S == Stage::Four) {
+        // compute the mean of dominance_checks_per_bucket
+        double mean_dominance_checks = 0.0;
+        for (size_t i = 0; i < dominance_checks_per_bucket.size(); ++i) {
+            mean_dominance_checks += dominance_checks_per_bucket[i];
+        }
+        auto step_calc = mean_dominance_checks / non_dominated_labels_per_bucket;
+        if (step_calc > 500) {
+
+            fmt::print("Increment redefinition counter\n");
+            if (redefine_counter % 5 == 0) {
+                redefine_counter = 0;
+                print_info("Step size is high, should increment step size - {}\n", step_calc);
+                fmt::print("Current step size bucket_interval: {}\n", bucket_interval);
+                redefine(bucket_interval + 20);
+                fmt::print("New step size bucket_interval: {}\n", bucket_interval);
+                reset_fixed_buckets();
+                fixed = false;
+            }
+            redefine_counter++;
+        }
 
+        // redefine_counter++;
+    }
     // Return the final list of merged labels after processing
     return merged_labels;
 }
@@ -786,8 +818,8 @@ inline bool BucketGraph::is_dominated(const Label *new_label, const Label *label
 #ifndef UNREACHABLE_DOMINANCE
     // Check visited jobs (bitmap comparison) for Stages 3, 4, and Enumerate
     if constexpr (S == Stage::Three || S == Stage::Four || S == Stage::Enumerate) {
-        // Iterate through the visited bitmap and ensure that the new label visits all jobs that the comparison label
-        // visits
+        // Iterate through the visited bitmap and ensure that the new label visits all jobs that the comparison
+        // label visits
         for (size_t i = 0; i < label->visited_bitmap.size(); ++i) {
             // If the comparison label visits a job that the new label does not, it is not dominated
             if ((label->visited_bitmap[i] & ~new_label->visited_bitmap[i]) != 0) { return false; }
diff --git a/bucket/BucketUtils.h b/bucket/BucketUtils.h
index e3c21d81..cac58107 100644
--- a/bucket/BucketUtils.h
+++ b/bucket/BucketUtils.h
@@ -64,80 +64,18 @@ void BucketGraph::add_arc(int from_bucket, int to_bucket, const std::vector<doub
  */
 template <Direction D>
 inline int BucketGraph::get_bucket_number(int job, const std::vector<double> &resource_values_vec) noexcept {
-    const auto &num_buckets_index = assign_buckets<D>(num_buckets_index_fw, num_buckets_index_bw);
-    const int   start_index       = num_buckets_index[job];
-
-    if constexpr (R_SIZE > 1) {
-        auto            &vrpjob = jobs[job];
-        std::vector<int> base_intervals(MAIN_RESOURCES);
-        std::vector<int> first_lb(MAIN_RESOURCES);
-        std::vector<int> first_ub(MAIN_RESOURCES);
-        std::vector<int> remainders(MAIN_RESOURCES);
-
-        // Compute base intervals, remainders, and bounds for each dimension
-        for (int r = 0; r < MAIN_RESOURCES; ++r) {
-            base_intervals[r] = (R_max[r] - R_min[r] + 1) / intervals[r].interval;
-            remainders[r]     = static_cast<int>(R_max[r] - R_min[r] + 1) % intervals[r].interval;
-            first_lb[r]       = vrpjob.lb[r]; // Lower bound for each dimension
-            first_ub[r]       = vrpjob.ub[r]; // Upper bound for each dimension
-        }
-
-        int bucket_index = 0;
-        int multiplier   = 1; // This will be used to compute the final bucket index across dimensions
-
-        if constexpr (D == Direction::Forward) {
-            // Iterate over each resource value and compute the bucket index for forward direction
-            for (int r = MAIN_RESOURCES - 1; r >= 0; --r) {
-                int resource_value_int = static_cast<int>(resource_values_vec[r]);
-                int bucket_dim_index;
-
-                // Handle uneven divisions by accounting for remainders
-                if (resource_value_int < first_lb[r] + base_intervals[r] * (remainders[r] > 0)) {
-                    bucket_dim_index = (resource_value_int - first_lb[r]) / base_intervals[r];
-                } else {
-                    bucket_dim_index = (resource_value_int - first_lb[r] - remainders[r]) / base_intervals[r];
-                }
-
-                bucket_index += bucket_dim_index * multiplier;
-                multiplier *=
-                    (R_max[r] - R_min[r] + 1) / base_intervals[r]; // Update the multiplier for the next dimension
-            }
-        } else if constexpr (D == Direction::Backward) {
-            // Iterate over each resource value and compute the bucket index for backward direction
-            for (int r = MAIN_RESOURCES - 1; r >= 0; --r) {
-                int resource_value_int = static_cast<int>(resource_values_vec[r]);
-                int bucket_dim_index;
-
-                // Handle uneven divisions by accounting for remainders
-                if (resource_value_int > first_ub[r] - base_intervals[r] * (remainders[r] > 0)) {
-                    bucket_dim_index = (first_ub[r] - resource_value_int) / base_intervals[r];
-                } else {
-                    bucket_dim_index = (first_ub[r] - resource_value_int - remainders[r]) / base_intervals[r];
-                }
-
-                bucket_index += bucket_dim_index * multiplier;
-                multiplier *=
-                    (R_max[r] - R_min[r] + 1) / base_intervals[r]; // Update the multiplier for the next dimension
-            }
-        }
-        return bucket_index + start_index; // Add the starting index to get the final bucket index
-    } else {
-        const auto &num_buckets_index = assign_buckets<D>(num_buckets_index_fw, num_buckets_index_bw);
-        const int   start_index       = num_buckets_index[job];
-        const int   base_interval     = (R_max[0] - R_min[0] + 1) / intervals[0].interval;
-        int         value             = static_cast<int>(resource_values_vec[0]);
-
-        const auto &vrpjob = jobs[job];
 
-        if constexpr (D == Direction::Forward) {
-            const int first_lb = vrpjob.lb[0];
-            return (value - first_lb) / base_interval + start_index;
-        } else if constexpr (D == Direction::Backward) {
-            const int first_ub = vrpjob.ub[0];
-            return (first_ub - value) / base_interval + start_index;
-        }
-        return -1; // If no bucket is found
+    std::vector<int> resource_vec_int(resource_values_vec.data(),
+                                      resource_values_vec.data() + resource_values_vec.size());
+    if constexpr (D == Direction::Forward) {
+        auto val = fw_job_interval_trees[job].query(resource_vec_int);
+        return val;
+    } else if constexpr (D == Direction::Backward) {
+        auto val = bw_job_interval_trees[job].query(resource_vec_int);
+        return val;
     }
+    std::throw_with_nested(std::runtime_error("BucketGraph::get_bucket_number: Invalid direction"));
+    return -1; // If no bucket is found
 }
 
 /**
@@ -159,67 +97,95 @@ void BucketGraph::define_buckets() {
     for (int r = 0; r < num_intervals; ++r) {
         total_ranges[r]   = R_max[r] - R_min[r] + 1;
         base_intervals[r] = total_ranges[r] / intervals[r].interval;
-        remainders[r]     = total_ranges[r] % intervals[r].interval; // Use std::fmod for floating-point modulo
+        remainders[r]     = total_ranges[r] % intervals[r].interval;
     }
+
     auto &buckets           = assign_buckets<D>(fw_buckets, bw_buckets);
     auto &num_buckets       = assign_buckets<D>(num_buckets_fw, num_buckets_bw);
     auto &num_buckets_index = assign_buckets<D>(num_buckets_index_fw, num_buckets_index_bw);
     num_buckets.resize(jobs.size());
     num_buckets_index.resize(jobs.size());
 
-    int cum_sum      = 0; // Keeps track of the global bucket index
-    int bucket_index = 0; // Keeps track of where to insert the next bucket
+    int cum_sum      = 0; // Tracks global bucket index
+    int bucket_index = 0;
+
+    // Helper function to check if the job can fit in a single bucket
+    auto fits_single_bucket = [&](const auto &job_total_ranges) {
+        for (int r = 0; r < num_intervals; ++r) {
+            if (job_total_ranges[r] > base_intervals[r]) { return false; }
+        }
+        return true;
+    };
 
     // Loop through each job to define its specific buckets
     for (const auto &VRPJob : jobs) {
         std::vector<int> job_total_ranges(num_intervals);
         for (int r = 0; r < num_intervals; ++r) { job_total_ranges[r] = VRPJob.ub[r] - VRPJob.lb[r]; }
-        // Check if the job's range fits within a single base interval for all dimensions
-        bool fits_single_bucket = true;
-        for (int r = 0; r < num_intervals; ++r) {
-            if (job_total_ranges[r] > base_intervals[r]) {
-                fits_single_bucket = false;
-                break;
-            }
-        }
 
-        if (fits_single_bucket) {
-            // Single bucket case for all resources
-            std::vector<int> lb(num_intervals), ub(num_intervals);
-            for (int r = 0; r < num_intervals; ++r) {
-                lb[r] = VRPJob.lb[r];
-                ub[r] = VRPJob.ub[r];
+        SplayTree job_tree;
+
+        if (fits_single_bucket(job_total_ranges)) {
+            // Single bucket case
+            std::vector<int> lb = VRPJob.lb, ub = VRPJob.ub;
+            buckets[bucket_index] = Bucket(VRPJob.id, lb, ub);
+
+            if constexpr (D == Direction::Forward) {
+                job_tree.insert(lb, ub, bucket_index);
+            } else {
+                job_tree.insert(lb, ub, bucket_index);
             }
-            buckets[bucket_index]        = Bucket(VRPJob.id, lb, ub);
+
             num_buckets[VRPJob.id]       = 1;
             num_buckets_index[VRPJob.id] = cum_sum;
             bucket_index++;
             cum_sum++;
+
         } else {
-            // Multiple bucket case, need to consider intervals for each resource dimension
+            // Multiple buckets case
             int              n_buckets = 0;
             std::vector<int> current_pos(num_intervals, 0);
 
-            // Nested loop over all intervals in each resource dimension
             while (true) {
                 std::vector<int> interval_start(num_intervals), interval_end(num_intervals);
 
+                // Calculate the start and end for each interval dimension
                 for (int r = 0; r < num_intervals; ++r) {
                     interval_start[r] = VRPJob.lb[r] + current_pos[r] * base_intervals[r];
                     interval_end[r]   = VRPJob.ub[r] - current_pos[r] * base_intervals[r];
 
+                    // Ensure intervals align with direction (forward/backward)
                     if constexpr (D == Direction::Forward) {
                         interval_end[r] = VRPJob.ub[r];
                     } else if constexpr (D == Direction::Backward) {
                         interval_start[r] = VRPJob.lb[r];
                     }
                 }
+
+                // Create a new bucket for this job
                 buckets[bucket_index] = Bucket(VRPJob.id, interval_start, interval_end);
+
+                if constexpr (D == Direction::Forward) {
+                    std::vector<int> interval_fw_end(interval_start.size());
+                    for (int r = 0; r < interval_start.size(); ++r) {
+                        interval_fw_end[r] = interval_start[r] + base_intervals[r];
+                    }
+
+                    job_tree.insert(interval_start, interval_fw_end, bucket_index);
+                } else {
+                    std::vector<int> interval_bw_end(interval_start.size());
+
+                    for (int r = 0; r < interval_start.size(); ++r) {
+                        interval_bw_end[r] = interval_end[r] - base_intervals[r];
+                    }
+
+                    job_tree.insert(interval_bw_end, interval_end, bucket_index);
+                }
+
                 bucket_index++;
                 n_buckets++;
                 cum_sum++;
 
-                // Update the position of intervals
+                // Check if we need to exit the loop (all intervals covered)
                 bool done = true;
                 for (int r = num_intervals - 1; r >= 0; --r) {
                     current_pos[r]++;
@@ -233,11 +199,18 @@ void BucketGraph::define_buckets() {
                 if (done) break;
             }
 
+            // Update job-specific bucket data
             num_buckets[VRPJob.id]       = n_buckets;
             num_buckets_index[VRPJob.id] = cum_sum - n_buckets;
         }
+        if constexpr (D == Direction::Forward) {
+            fw_job_interval_trees[VRPJob.id] = job_tree;
+        } else {
+            bw_job_interval_trees[VRPJob.id] = job_tree;
+        }
     }
 
+    // Update global bucket sizes based on direction
     if constexpr (D == Direction::Forward) {
         fw_buckets_size = cum_sum;
     } else {
@@ -357,7 +330,8 @@ void BucketGraph::generate_arcs() {
                     std::lock_guard<std::mutex> lock(buckets_mutex); // Lock the mutex to ensure thread safety
                     add_arc<D>(from_bucket, to_bucket, local_res_inc, local_cost_inc); // Add the arc globally
                     buckets[from_bucket].add_bucket_arc(from_bucket, to_bucket, local_res_inc, local_cost_inc,
-                                                        D == Direction::Forward, false); // Add the arc to the bucket
+                                                        D == Direction::Forward,
+                                                        false); // Add the arc to the bucket
                 }
             }
         }
diff --git a/examples/VRPTW.cpp b/examples/VRPTW.cpp
index d356e4f2..4a0f74ac 100644
--- a/examples/VRPTW.cpp
+++ b/examples/VRPTW.cpp
@@ -206,8 +206,8 @@ int main(int argc, char *argv[]) {
         jobs[k].set_location(instance.x_coord[k], instance.y_coord[k]);
         jobs[k].lb.push_back(start_time);
         jobs[k].ub.push_back(end_time);
-        jobs[k].lb.push_back(0);
-        jobs[k].ub.push_back(instance.q);
+        //jobs[k].lb.push_back(0);
+        //jobs[k].ub.push_back(instance.q);
         jobs[k].consumption.push_back(duration);
         jobs[k].consumption.push_back(demand);
     }
diff --git a/examples/VRPTW.h b/examples/VRPTW.h
index a395b632..a310bba5 100644
--- a/examples/VRPTW.h
+++ b/examples/VRPTW.h
@@ -831,7 +831,7 @@ class VRProblem {
                     ////////////////////////////////////////////////////
                     // Handle non-violated cuts in a single pass
                     ////////////////////////////////////////////////////
-                    bool cleared     = false;
+                    bool cleared        = false;
                     auto n_cuts_removed = 0;
                     // Iterate over the constraints in reverse order to remove non-violated cuts
                     for (int i = SRCconstraints.size() - 1; i >= 0; --i) {
diff --git a/external/hgs_vrptw/Individual.cpp b/external/hgs_vrptw/Individual.cpp
index e7f8e3d7..c7181db8 100644
--- a/external/hgs_vrptw/Individual.cpp
+++ b/external/hgs_vrptw/Individual.cpp
@@ -303,8 +303,8 @@ Individual::Individual(Params *params, bool rcws, std::vector<std::vector<int>>
                             }
                             chromR[nextEmptyRoute].push_back(tournamentSavings[i].c1);
                             chromR[nextEmptyRoute].push_back(tournamentSavings[i].c2);
-                            load[nextEmptyRoute] +=
-                                params->cli[tournamentSavings[i].c1].demand + params->cli[tournamentSavings[i].c2].demand;
+                            load[nextEmptyRoute] += params->cli[tournamentSavings[i].c1].demand +
+                                                    params->cli[tournamentSavings[i].c2].demand;
                             inRoute[tournamentSavings[i].c1] = true;
                             inRoute[tournamentSavings[i].c2] = true;
                             while (nextEmptyRoute < nbVehicles && !chromR[nextEmptyRoute].empty()) nextEmptyRoute++;
diff --git a/external/hgs_vrptw/Individual.h b/external/hgs_vrptw/Individual.h
index a901df5a..5157c309 100644
--- a/external/hgs_vrptw/Individual.h
+++ b/external/hgs_vrptw/Individual.h
@@ -30,14 +30,13 @@ SOFTWARE.*/
 
 #include "Params.h"
 
-struct EvalIndiv
-{
-	double penalizedCost = 0.;		// Penalized cost of the solution
-	int nbRoutes = 0;				// Number of routes
-	double distance = 0.;			// Total distance
-	double capacityExcess = 0.;		// Sum of excess load in all routes
-	double durationExcess = 0.;		// Sum of excess duration in all routes
-	bool isFeasible = false;		// Feasibility status of the individual
+struct EvalIndiv {
+    double penalizedCost  = 0.;    // Penalized cost of the solution
+    int    nbRoutes       = 0;     // Number of routes
+    double distance       = 0.;    // Total distance
+    double capacityExcess = 0.;    // Sum of excess load in all routes
+    double durationExcess = 0.;    // Sum of excess duration in all routes
+    bool   isFeasible     = false; // Feasibility status of the individual
 };
 
 // Object to store all relevant information that may be needed to calculate some cost corresponding to a solution
diff --git a/external/hgs_vrptw/Population.cpp b/external/hgs_vrptw/Population.cpp
index 73cf7a81..7984c357 100644
--- a/external/hgs_vrptw/Population.cpp
+++ b/external/hgs_vrptw/Population.cpp
@@ -643,4 +643,4 @@ void Population::mineElite() {
         mdmEliteUpdated = false;
         mdmNextPattern  = 0;
     }
-}
\ No newline at end of file
+}
diff --git a/include/Definitions.h b/include/Definitions.h
index c0d97dbc..51efb190 100644
--- a/include/Definitions.h
+++ b/include/Definitions.h
@@ -1096,3 +1096,353 @@ inline ModelData extractModelData(GRBModel &model) {
 
     return data;
 }
+
+class IntervalTreeNode {
+public:
+    int               low, high;    // Interval [low, high]
+    int               bucket_index; // Associated bucket index
+    int               job_id;       // Job ID associated with this interval
+    int               max;          // Max endpoint in the subtree
+    IntervalTreeNode *left;
+    IntervalTreeNode *right;
+
+    // Constructor
+    IntervalTreeNode(int l, int h, int idx, int job)
+        : low(l), high(h), bucket_index(idx), job_id(job), max(h), left(nullptr), right(nullptr) {}
+};
+
+class IntervalTree {
+public:
+    IntervalTreeNode *root;
+
+    IntervalTree() : root(nullptr) {}
+
+    // Insert interval [low, high] with bucket_index and job_id
+    IntervalTreeNode *insert(IntervalTreeNode *node, int low, int high, int bucket_index, int job_id) {
+        if (!node) { return new IntervalTreeNode(low, high, bucket_index, job_id); }
+
+        // Inserting based on the low value
+        if (low < node->low) {
+            node->left = insert(node->left, low, high, bucket_index, job_id);
+        } else if (low > node->low) {
+            node->right = insert(node->right, low, high, bucket_index, job_id);
+        } else {
+            // Handle the case where two intervals have the same 'low' value.
+            // We might choose to insert this to the right if 'high' values differ.
+            if (high > node->high) {
+                node->right = insert(node->right, low, high, bucket_index, job_id);
+            } else {
+                node->left = insert(node->left, low, high, bucket_index, job_id);
+            }
+        }
+
+        // Update the max value of this node
+        node->max = std::max(node->max, high);
+
+        return node;
+    }
+
+    // Query to find the bucket index for a given resource value
+    int query(IntervalTreeNode *node, int value) {
+        if (!node) {
+            std::cout << "Query failed: no matching node for value " << value << std::endl;
+            // print the tree
+            for (int i = 0; i < 10; ++i) {
+                std::cout << "Node " << i << ": ";
+                if (node) {
+                    std::cout << "[" << node->low << ", " << node->high << "] with max " << node->max << std::endl;
+                } else {
+                    std::cout << "NULL" << std::endl;
+                }
+            }
+            return -1;
+        }
+
+        // Check if the value lies within or is exactly on the boundary of this interval
+        if (value >= node->low && value <= node->high) { return node->bucket_index; }
+
+        // Move to left subtree if the value is less than or equal to the max of the left subtree
+        if (node->left && value <= node->left->max) { return query(node->left, value); }
+
+        // Otherwise, move to the right subtree
+        return query(node->right, value);
+    }
+
+    // Helper to insert interval into the tree
+    void insert(int low, int high, int bucket_index, int job_id) {
+        root = insert(root, low, high, bucket_index, job_id);
+    }
+
+    // Helper to query the tree
+    int query(int value) { return query(root, value); }
+};
+
+/**
+ * @class SplayTree
+ * @brief A class representing a Splay Tree, which is a self-adjusting binary search tree.
+ *
+ * The Splay Tree supports efficient insertion, deletion, and search operations by performing
+ * splay operations that move accessed nodes closer to the root, thereby improving access times
+ * for frequently accessed nodes.
+ *
+ * The tree nodes store intervals [low, high] and a bucket index associated with each interval.
+ *
+ * @note This implementation assumes that the TreeNode class is defined elsewhere with the
+ *       necessary members: low, high, bucket_index, left, right, and parent.
+
+ */
+class TreeNode {
+public:
+    std::vector<int> low;          // Lower bounds for each dimension
+    std::vector<int> high;         // Upper bounds for each dimension
+    int              bucket_index; // Bucket index for this node
+    TreeNode        *left;
+    TreeNode        *right;
+    TreeNode        *parent;
+
+    TreeNode(const std::vector<int> &low, const std::vector<int> &high, int bucket_index)
+        : low(low), high(high), bucket_index(bucket_index), left(nullptr), right(nullptr), parent(nullptr) {}
+
+    bool contains(const std::vector<int> &point) const {
+        for (size_t i = 0; i < low.size(); ++i) {
+            if (point[i] < low[i] || point[i] > high[i]) { return false; }
+        }
+        return true;
+    }
+
+    bool is_less_than(const std::vector<int> &point) const {
+        for (size_t i = 0; i < low.size(); ++i) {
+            if (high[i] < point[i]) {
+                return true;
+            } else if (low[i] > point[i]) {
+                return false;
+            }
+        }
+        return false; // This case shouldn't be reached if comparing proper intervals
+    }
+};
+
+class SplayTree {
+    TreeNode *root;
+
+    void zig(TreeNode *x) {
+        TreeNode *p = x->parent;
+        TreeNode *B = (p->left == x) ? x->right : x->left;
+
+        x->parent = p->parent;
+        p->parent = x;
+
+        if (p->left == x) {
+            x->right = p;
+            p->left  = B;
+        } else {
+            x->left  = p;
+            p->right = B;
+        }
+
+        if (B != nullptr) B->parent = p;
+    }
+
+    void zig_zig(TreeNode *x) {
+        TreeNode *p = x->parent;
+        TreeNode *g = p->parent;
+        if (p->left == x) {
+            TreeNode *B = x->right;
+            TreeNode *C = p->right;
+
+            x->parent = g->parent;
+            x->right  = p;
+
+            p->parent = x;
+            p->left   = B;
+            p->right  = g;
+
+            g->parent = p;
+            g->left   = C;
+
+            if (x->parent != nullptr) {
+                if (x->parent->left == g)
+                    x->parent->left = x;
+                else
+                    x->parent->right = x;
+            }
+
+            if (B != nullptr) B->parent = p;
+            if (C != nullptr) C->parent = g;
+        } else {
+            TreeNode *B = p->left;
+            TreeNode *C = x->left;
+
+            x->parent = g->parent;
+            x->left   = p;
+
+            p->parent = x;
+            p->left   = g;
+            p->right  = C;
+
+            g->parent = p;
+            g->right  = B;
+
+            if (x->parent != nullptr) {
+                if (x->parent->left == g)
+                    x->parent->left = x;
+                else
+                    x->parent->right = x;
+            }
+
+            if (B != nullptr) B->parent = g;
+            if (C != nullptr) C->parent = p;
+        }
+    }
+
+    void zig_zag(TreeNode *x) {
+        TreeNode *p = x->parent;
+        TreeNode *g = p->parent;
+
+        if (p->right == x) {
+            TreeNode *B = x->left;
+            TreeNode *C = x->right;
+
+            x->parent = g->parent;
+            x->left   = p;
+            x->right  = g;
+
+            p->parent = x;
+            p->right  = B;
+
+            g->parent = x;
+            g->left   = C;
+
+            if (x->parent != nullptr) {
+                if (x->parent->left == g)
+                    x->parent->left = x;
+                else
+                    x->parent->right = x;
+            }
+
+            if (B != nullptr) B->parent = p;
+            if (C != nullptr) C->parent = g;
+        } else {
+            TreeNode *B = x->left;
+            TreeNode *C = x->right;
+
+            x->parent = g->parent;
+            x->left   = g;
+            x->right  = p;
+
+            p->parent = x;
+            p->left   = C;
+
+            g->parent = x;
+            g->right  = B;
+
+            if (x->parent != nullptr) {
+                if (x->parent->left == g)
+                    x->parent->left = x;
+                else
+                    x->parent->right = x;
+            }
+
+            if (B != nullptr) B->parent = g;
+            if (C != nullptr) C->parent = p;
+        }
+    }
+
+    void splay(TreeNode *x) {
+        while (x->parent != nullptr) {
+            TreeNode *p = x->parent;
+            TreeNode *g = p->parent;
+
+            if (g == nullptr) {
+                zig(x);
+            } else if (g->left == p && p->left == x) {
+                zig_zig(x);
+            } else if (g->right == p && p->right == x) {
+                zig_zig(x);
+            } else {
+                zig_zag(x);
+            }
+        }
+        root = x;
+    }
+
+public:
+    SplayTree() : root(nullptr) {}
+
+    TreeNode *find(const std::vector<int> &point) {
+        TreeNode *curr = root;
+        TreeNode *prev = nullptr;
+
+        while (curr != nullptr) {
+            prev = curr;
+            if (curr->contains(point)) {
+                splay(curr);
+                return curr;
+            } else if (curr->is_less_than(point)) {
+                curr = curr->right;
+            } else {
+                curr = curr->left;
+            }
+        }
+
+        if (prev != nullptr) splay(prev);
+        return nullptr;
+    }
+
+    int query(const std::vector<int> &point) {
+        TreeNode *node = find(point);
+        if (node != nullptr) return node->bucket_index;
+        return -1;
+    }
+
+    // Insert a new multidimensional interval
+    void insert(const std::vector<int> &low, const std::vector<int> &high, int bucket_index) {
+
+        if (root == nullptr) {
+            root = new TreeNode(low, high, bucket_index);
+            return;
+        }
+
+        TreeNode *curr = root;
+        while (curr != nullptr) {
+            if (low < curr->low) {
+                if (curr->left == nullptr) {
+                    TreeNode *newNode = new TreeNode(low, high, bucket_index);
+                    curr->left        = newNode;
+                    newNode->parent   = curr;
+                    splay(newNode);
+                    return;
+                } else {
+                    curr = curr->left;
+                }
+            } else if (low > curr->low) {
+                if (curr->right == nullptr) {
+                    TreeNode *newNode = new TreeNode(low, high, bucket_index);
+                    curr->right       = newNode;
+                    newNode->parent   = curr;
+                    splay(newNode);
+                    return;
+                } else {
+                    curr = curr->right;
+                }
+            } else {
+                splay(curr);
+                return; // Duplicate interval
+            }
+        }
+    }
+
+    void inOrderPrint(TreeNode *node) {
+        if (node == nullptr) return;
+        inOrderPrint(node->left);
+        std::cout << "Bucket " << node->bucket_index << ": [";
+        for (size_t i = 0; i < node->low.size(); ++i) {
+            std::cout << "[" << node->low[i] << ", " << node->high[i] << "]";
+            if (i != node->low.size() - 1) std::cout << ", ";
+        }
+        std::cout << "]\n";
+        inOrderPrint(node->right);
+    }
+
+    void print() { inOrderPrint(root); }
+};
diff --git a/include/HGS.h b/include/HGS.h
index 111b1ee3..aaf5c67b 100644
--- a/include/HGS.h
+++ b/include/HGS.h
@@ -37,16 +37,16 @@ class HGS {
                   << " VEHICLES" << std::endl;
         std::cout << "----- BUILDING INITIAL POPULATION" << std::endl;
         Population population(&params, &split, &localSearch);
-
-        
-        // Genetic algorithm
-        std::cout << "----- STARTING GENETIC ALGORITHM" << std::endl;
-        Genetic solver(&params, &split, &population, &localSearch);
-        solver.run(config.nbIter, config.timeLimit);
-        std::cout << "----- GENETIC ALGORITHM FINISHED, TIME SPENT: " << params.getTimeElapsedSeconds() << std::endl;
-
+        /*
+                // Genetic algorithm
+                std::cout << "----- STARTING GENETIC ALGORITHM" << std::endl;
+                Genetic solver(&params, &split, &population, &localSearch);
+                solver.run(config.nbIter, config.timeLimit);
+                std::cout << "----- GENETIC ALGORITHM FINISHED, TIME SPENT: " << params.getTimeElapsedSeconds() <<
+           std::endl;
+                */
         auto sol = population.extractFeasibleRoutes();
-        
+
         // Return 0 if the program execution was successfull
         return sol;
     }
diff --git a/src/BucketGraph.cpp b/src/BucketGraph.cpp
index 8c8dbc6c..b930635f 100644
--- a/src/BucketGraph.cpp
+++ b/src/BucketGraph.cpp
@@ -357,7 +357,7 @@ void BucketGraph::calculate_neighborhoods(size_t num_closest) {
     size_t num_jobs = jobs.size();
 
     // Initialize the neighborhood bitmaps as vectors of uint64_t for forward and backward neighborhoods
-    neighborhoods_bitmap.resize(num_jobs);                           // Forward neighborhood
+    neighborhoods_bitmap.resize(num_jobs); // Forward neighborhood
 
     for (size_t i = 0; i < num_jobs; ++i) {
         std::vector<std::pair<double, size_t>> forward_distances; // Distances for forward neighbors
@@ -627,9 +627,15 @@ void BucketGraph::common_initialization() {
     fw_c_bar.clear();
     bw_c_bar.clear();
 
+    // print fw_buckets size
+    // fmt::print("fw_buckets size: {}\n", fw_buckets.size());
+    // fmt::print("fw_buckets size: {}\n", bw_buckets_size);
+    dominance_checks_per_bucket.assign(fw_buckets_size + 1, 0);
+    non_dominated_labels_per_bucket = 0;
+
     // Resize cost vectors to match the number of buckets
-    fw_c_bar.resize(fw_buckets.size(), std::numeric_limits<double>::infinity());
-    bw_c_bar.resize(bw_buckets.size(), std::numeric_limits<double>::infinity());
+    fw_c_bar.resize(fw_buckets_size, std::numeric_limits<double>::infinity());
+    bw_c_bar.resize(bw_buckets_size, std::numeric_limits<double>::infinity());
 
     auto &num_buckets      = assign_buckets<Direction::Forward>(num_buckets_fw, num_buckets_bw);
     auto &num_bucket_index = assign_buckets<Direction::Forward>(num_buckets_index_fw, num_buckets_index_bw);
@@ -647,7 +653,7 @@ void BucketGraph::common_initialization() {
     }
 
     // Clear forward and backward buckets
-    for (auto b = 0; b < fw_buckets_size; b++) {
+    for (auto b = 0; b < fw_buckets.size(); b++) {
         fw_buckets[b].clear();
         bw_buckets[b].clear();
     }
@@ -1129,4 +1135,3 @@ Label *BucketGraph::compute_mono_label(const Label *L) {
 
     return new_label;
 }
-