Broken Pairs Exchange & Subpopulations (#72)

* Introduce Subpopulation struct

* New statistics members for subpopulation diversity

* Zero diversity if no other population neighbors; fix accumulate

* Add broken pairs exchange

* Destroy consecutive broken pairs

* Only destroy/repair worst parent

* Streamline config, rename pop to subpop, refactor BPX

* Change nbIterNoImprove counter based on population current best

* Address small review comments

* Refactor subpopulations

* Make lists for crossover operators

* Fix indexing bug

* Removing comments, renaming variables

* Fix export indiv in analysis

* Re-add fitness calculation conditional

* Remove weird line, remove diversityWeight, change res.get_run_time

* Bugfix population nbIterNoImprove mechanism

* Type alias SubPopulation, change Member to IndividualWrapper

* Fix parent comparison in parent selection

* Refactor Individual::operator==

* Bugfix type and indexing BPX

* Replace Inidividual pointers with unique_ptr

* Rename getCurrentBest to getCurrentBestFeasibleCost

* Simplify updateBiasedFitness and population generation;

* Simplify fitness calculation

* Refactor getBinaryTournament

* Replace size check with feasible.empty

* Address review comments

* Rename subs to wrapper

analysis.py

@@ -77,9 +77,16 @@ def solve(loc: str, seed: int, **kwargs):
         ls.add_route_operator(op)
 
     algo = hgspy.GeneticAlgorithm(params, rng, pop, ls)
-    algo.add_crossover_operator(hgspy.crossover.alternating_exchange)
-    algo.add_crossover_operator(hgspy.crossover.ordered_exchange)
-    algo.add_crossover_operator(hgspy.crossover.selective_route_exchange)
+
+    crossover_ops = [
+        hgspy.crossover.alternating_exchange,
+        hgspy.crossover.broken_pairs_exchange,
+        hgspy.crossover.ordered_exchange,
+        hgspy.crossover.selective_route_exchange,
+    ]
+
+    for op in crossover_ops:
+        algo.add_crossover_operator(op)
 
     if "phase" in kwargs and kwargs["phase"]:
         t_lim = tools.static_time_limit(tools.name2size(loc), kwargs["phase"])
@@ -138,7 +145,7 @@ def make_path(subdir, extension):
     sol_path = make_path(_SOLS_DIR, "txt")
     best = res.get_best_found()
     stats = res.get_statistics()
-    best.export_cvrplib_format(sol_path, sum(stats.run_times()))
+    best.export_cvrplib_format(sol_path, res.get_run_time())
 
     # Save statistics
     stats_path = make_path(_STATS_DIR, "csv")

benchmark.py

@@ -63,9 +63,16 @@ def solve(loc: str, seed: int, **kwargs):
         ls.add_route_operator(op)
 
     algo = hgspy.GeneticAlgorithm(params, rng, pop, ls)
-    algo.add_crossover_operator(hgspy.crossover.alternating_exchange)
-    algo.add_crossover_operator(hgspy.crossover.ordered_exchange)
-    algo.add_crossover_operator(hgspy.crossover.selective_route_exchange)
+
+    crossover_ops = [
+        hgspy.crossover.alternating_exchange,
+        hgspy.crossover.broken_pairs_exchange,
+        hgspy.crossover.ordered_exchange,
+        hgspy.crossover.selective_route_exchange,
+    ]
+
+    for op in crossover_ops:
+        algo.add_crossover_operator(op)
 
     if "phase" in kwargs and kwargs["phase"]:
         t_lim = tools.static_time_limit(tools.name2size(loc), kwargs["phase"])

hgs_vrptw/include/CommandLine.h

@@ -80,8 +80,6 @@ class CommandLine
             else if (std::string(argv[i]) == "-intensificationProbability")
                 config.intensificationProbability
                     = static_cast<size_t>(atoi(argv[i + 1]));
-            else if (std::string(argv[i]) == "-diversityWeight")
-                config.diversityWeight = atof(argv[i + 1]);
             else if (std::string(argv[i])
                      == "-circleSectorOverlapToleranceDegrees")
             {

hgs_vrptw/include/Config.h

@@ -35,10 +35,6 @@ struct Config
 
     size_t selectProbability = 90;  // offspring selection probability
 
-    // Weight for diversity criterion. If 0, weight is set to 1 - nbElite /
-    // populationSize
-    double diversityWeight = 0.0;
-
     int nbVeh = INT_MAX;  // Number of vehicles
 
     // Granular search parameter, limits the number of moves in the RI local
@@ -60,6 +56,9 @@ struct Config
     // that even small circle sectors have 'overlap'
     int minCircleSectorSize = static_cast<int>(15 / 360. * 65536);
 
+    // Percentage of customers to remove in brokenPairsExchange
+    size_t destroyPct = 20;
+
     explicit Config(int seed = 0,
                     size_t nbIter = 20'000,
                     int timeLimit = INT_MAX,
@@ -77,14 +76,14 @@ struct Config
                     size_t repairProbability = 50,
                     size_t repairBooster = 10,
                     size_t selectProbability = 90,
-                    double diversityWeight = 0.,
                     int nbVeh = INT_MAX,
                     size_t nbGranular = 40,
                     int weightWaitTime = 2,
                     int weightTimeWarp = 10,
                     size_t intensificationProbability = 25,
                     int circleSectorOverlapToleranceDegrees = 0,
-                    int minCircleSectorSizeDegrees = 15)
+                    int minCircleSectorSizeDegrees = 15,
+                    size_t destroyPct = 20)
         : seed(seed),
           nbIter(nbIter),
           timeLimit(timeLimit),
@@ -102,12 +101,12 @@ struct Config
           repairProbability(repairProbability),
           repairBooster(repairBooster),
           selectProbability(selectProbability),
-          diversityWeight(diversityWeight),
           nbVeh(nbVeh),
           nbGranular(nbGranular),
           weightWaitTime(weightWaitTime),
           weightTimeWarp(weightTimeWarp),
-          intensificationProbability(intensificationProbability)
+          intensificationProbability(intensificationProbability),
+          destroyPct(destroyPct)
     {
         auto const overlap = circleSectorOverlapToleranceDegrees / 360. * 65536;
         circleSectorOverlapTolerance = static_cast<int>(overlap);

hgs_vrptw/include/Individual.h

@@ -126,6 +126,11 @@ class Individual
         return cost() < other.cost();
     }
 
+    bool operator==(Individual const &other) const
+    {
+        return cost() == other.cost() && routes_ == other.routes_;
+    }
+
     Individual(Params const *params, XorShift128 *rng);  // random individual
 
     Individual(Params const *params, Tour tour);

hgs_vrptw/include/Population.h

@@ -6,6 +6,7 @@
 #include "Statistics.h"
 #include "XorShift128.h"
 
+#include <memory>
 #include <vector>
 
 // Class representing the population of a genetic algorithm with do binary
@@ -14,25 +15,33 @@ class Population
 {
     friend class Statistics;  // used to collect population statistics
 
+    struct IndividualWrapper
+    {
+        std::unique_ptr<Individual> indiv;
+        size_t fitness;
+    };
+
+    using SubPopulation = std::vector<IndividualWrapper>;
     using Parents = std::pair<Individual const *, Individual const *>;
 
     Params &params;    // Problem parameters
     XorShift128 &rng;  // Random number generator
 
-    std::vector<Individual *> population;  // Population ordered asc. by cost
-    std::vector<double> fitness;           // Population fitness
-    Individual bestSol;                    // best observed solution
+    SubPopulation feasible;    // Sub-population ordered asc. by cost
+    SubPopulation infeasible;  // Sub-population ordered asc. by cost
+
+    Individual bestSol;
 
-    // Evaluates the biased fitness of all individuals in the population
-    void updateBiasedFitness();
+    // Evaluates the biased fitness of all individuals in the sub-population
+    void updateBiasedFitness(SubPopulation &subPop);
 
-    // Removes a duplicate individual from the population if there exists one.
-    // If there are multiple duplicate individuals, then the one with the lowest
-    // index in `population` is removed first.
-    bool removeDuplicate();
+    // Removes a duplicate individual from the sub-population if there exists
+    // one. If there are multiple duplicate individuals, then the one with the
+    // lowest index in the sub-population is removed first.
+    bool removeDuplicate(SubPopulation &subPop);
 
     // Removes the worst individual in terms of biased fitness
-    void removeWorstBiasedFitness();
+    void removeWorstBiasedFitness(SubPopulation &subPop);
 
     // Generates a population of passed-in size
     void generatePopulation(size_t popSize);
@@ -55,11 +64,11 @@ class Population
      */
     void reorder()
     {
-        std::sort(population.begin(),
-                  population.end(),
-                  [](auto const &indiv1, auto const &indiv2) {
-                      return indiv1->cost() < indiv2->cost();
-                  });
+        auto const op = [](auto const &wrapper1, auto const &wrapper2) {
+            return wrapper1.indiv->cost() < wrapper2.indiv->cost();
+        };
+        std::sort(feasible.begin(), feasible.end(), op);
+        std::sort(infeasible.begin(), infeasible.end(), op);
     }
 
     // Selects two (if possible non-identical) parents by binary tournament
@@ -70,9 +79,16 @@ class Population
      */
     [[nodiscard]] Individual const &getBestFound() const { return bestSol; }
 
-    Population(Params &params, XorShift128 &rng);
+    /**
+     * Returns the current best objective value in the feasible sub-population
+     * or ``INT_MAX`` if no feasible solution exists.
+     */
+    [[nodiscard]] size_t getCurrentBestFeasibleCost() const
+    {
+        return !feasible.empty() ? feasible[0].indiv->cost() : INT_MAX;
+    }
 
-    ~Population();
+    Population(Params &params, XorShift128 &rng);
 };
 
 #endif

hgs_vrptw/include/Statistics.h

@@ -21,13 +21,14 @@ class Statistics
     std::vector<double> iterTimes_;
     std::vector<size_t> popSizes_;
     std::vector<size_t> numFeasiblePop_;
-    std::vector<double> popDiversity_;
-    std::vector<size_t> penaltiesCapacity_;
-    std::vector<size_t> penaltiesTimeWarp_;
+    std::vector<double> feasDiversity_;
     std::vector<size_t> feasBest_;
     std::vector<size_t> feasAverage_;
+    std::vector<double> infeasDiversity_;
     std::vector<size_t> infeasBest_;
     std::vector<size_t> infeasAverage_;
+    std::vector<size_t> penaltiesCapacity_;
+    std::vector<size_t> penaltiesTimeWarp_;
 
     timedDatapoints incumbents_;
 
@@ -85,31 +86,15 @@ class Statistics
     }
 
     /**
-     * Returns a vector of the average population diversity, one element per
-     * iteration. The average diversity is computed as the average broken pairs
-     * distance for each individual in the population, compared to its
-     * neighbours (the neighbourhood size is controlled by the ``nbClose``
-     * setting).
+     * Returns a vector of the average feasible sub-population diversity, one
+     * element per iteration. The average diversity is computed as the average
+     * broken pairs distance for each individual in the sub-population, compared
+     * to its neighbours (the neighbourhood size is controlled by the
+     * ``nbClose`` setting).
      */
-    [[nodiscard]] std::vector<double> const &popDiversity() const
+    [[nodiscard]] std::vector<double> const &feasDiversity() const
     {
-        return popDiversity_;
-    }
-
-    /**
-     * Returns a vector of capacity penalties, one element per iteration.
-     */
-    [[nodiscard]] std::vector<size_t> const &penaltiesCapacity() const
-    {
-        return penaltiesCapacity_;
-    }
-
-    /**
-     * Returns a vector of time warp penalties, one element per iteration.
-     */
-    [[nodiscard]] std::vector<size_t> const &penaltiesTimeWarp() const
-    {
-        return penaltiesTimeWarp_;
+        return feasDiversity_;
     }
 
     /**
@@ -132,6 +117,18 @@ class Statistics
         return feasAverage_;
     }
 
+    /**
+     * Returns a vector of the average infeasible sub-population diversity, one
+     * element per iteration. The average diversity is computed as the average
+     * broken pairs distance for each individual in the sub-population, compared
+     * to its neighbours (the neighbourhood size is controlled by the
+     * ``nbClose`` setting).
+     */
+    [[nodiscard]] std::vector<double> const &infeasDiversity() const
+    {
+        return infeasDiversity_;
+    }
+
     /**
      * Returns a vector of the best objective value of infeasible individuals,
      * one element per iteration. If there are no infeasible individuals, then
@@ -152,6 +149,22 @@ class Statistics
         return infeasAverage_;
     }
 
+    /**
+     * Returns a vector of capacity penalties, one element per iteration.
+     */
+    [[nodiscard]] std::vector<size_t> const &penaltiesCapacity() const
+    {
+        return penaltiesCapacity_;
+    }
+
+    /**
+     * Returns a vector of time warp penalties, one element per iteration.
+     */
+    [[nodiscard]] std::vector<size_t> const &penaltiesTimeWarp() const
+    {
+        return penaltiesTimeWarp_;
+    }
+
     /**
      * Returns a vector of (runtime, objective)-pairs, one for each time
      * a new, feasible best heuristic solution has been found.

hgs_vrptw/include/crossover.h

@@ -57,4 +57,15 @@ Individual alternatingExchange(
     Params const &params,
     XorShift128 &rng);
 
+/**
+ * Performs one Broken Pair Crossover of the given parents. A client is removed
+ * from the worst parents if its successor is not identical to the client's
+ * sucessor in the other parent. Removed clients are greedily re-inserted in the
+ * solution.
+ */
+Individual brokenPairsExchange(
+    std::pair<Individual const *, Individual const *> const &parents,
+    Params const &params,
+    XorShift128 &rng);
+
 #endif  // CROSSOVER_H

hgs_vrptw/src/CMakeLists.txt

@@ -4,6 +4,7 @@ set(geneticSources
         GeneticAlgorithm.cpp
         crossover/crossover.cpp
         crossover/alternatingExchange.cpp
+        crossover/brokenPairsExchange.cpp
         crossover/orderedExchange.cpp
         crossover/selectiveRouteExchange.cpp)
 

hgs_vrptw/src/GeneticAlgorithm.cpp

@@ -25,28 +25,30 @@ Result GeneticAlgorithm::run(StoppingCriterion &stop)
     Statistics stats;
 
     size_t iter = 0;
-    size_t nbIterNonProd = 1;
+    size_t nbIterNoImprove = 1;
 
     auto start = clock::now();
     while (not stop())
     {
         iter++;
 
-        if (nbIterNonProd == params.config.nbIter)  // restart population after
-        {                                           // this number of useless
-            population.restart();                   // iterations
-            nbIterNonProd = 1;
+        if (nbIterNoImprove == params.config.nbIter)  // restart population
+        {                                             // after this number of
+            population.restart();                     // non-improving iters
+            nbIterNoImprove = 1;
         }
 
-        auto const currBest = population.getBestFound().cost();
+        auto const currBest = population.getCurrentBestFeasibleCost();
 
         auto offspring = crossover();
         educate(offspring);
 
-        if (currBest > population.getBestFound().cost())  // has new best!
-            nbIterNonProd = 1;
+        auto const newBest = population.getCurrentBestFeasibleCost();
+
+        if (currBest > newBest)  // has new best!
+            nbIterNoImprove = 1;
         else
-            nbIterNonProd++;
+            nbIterNoImprove++;
 
         // Diversification and penalty management
         if (iter % params.config.nbPenaltyManagement == 0)

hgs_vrptw/src/Individual.cpp

@@ -253,6 +253,9 @@ void Individual::brokenPairsDistance(Individual *other)
 
 double Individual::avgBrokenPairsDistanceClosest() const
 {
+    if (indivsPerProximity.empty())
+        return 0;
+
     size_t maxSize
         = std::min(params->config.nbClose, indivsPerProximity.size());