Adding discrete seeds to TSP solvers

docs/solvers.rst

@@ -111,13 +111,18 @@ Notice this local optimum may be different for distinct perturbation schemes and
 
     max_processing_time {None}
         Maximum processing time in seconds. If not provided, the method stops
-        only when a local minimum is obtained
+        only when a local minimum is obtained or if max_iterations is set.
 
-    log_file
+    rng {None}
+        A python random number generator which can be seeded to ensure 
+        consistent results given a call to this function. If not provided
+        the initial random number generator is used.
+
+    log_file {None}
         If not `None`, creates a log file with details about the whole
         execution
 
-    verbose
+    verbose {False}
         If true, prints algorithm status every iteration
 
 
@@ -138,6 +143,8 @@ An implementation of the `Simulated Annealing <https://en.wikipedia.org/wiki/Sim
       perturbation_scheme: str = "two_opt",
       alpha: float = 0.9,
       max_processing_time: Optional[float] = None,
+      max_iterations: Optional[int] = None,
+      rng: Optional[Random] = None,
       log_file: Optional[str] = None,
       verbose: bool = False,
   )
@@ -165,7 +172,13 @@ An implementation of the `Simulated Annealing <https://en.wikipedia.org/wiki/Sim
 
     max_processing_time {None}
         Maximum processing time in seconds. If not provided, the method stops
-        only when there were 3 temperature cycles with no improvement.
+        only when there were 3 temperature cycles with no improvement or if
+        max_iterations is set.
+
+    rng {None}
+        A python random number generator which can be seeded to ensure 
+        consistent results given a call to this function. If not provided
+        the initial random number generator is used.
 
     log_file {None}
         If not `None`, creates a log file with details about the whole
@@ -209,11 +222,11 @@ A basic Lin and Kernighan implementation is provided. It can be said that the qu
     x0
         Initial permutation. If not provided, it starts with a random path.
 
-    log_file
+    log_file {None}
         If not `None`, creates a log file with details about the whole
         execution.
 
-    verbose
+    verbose {False}
         If true, prints algorithm status every iteration.
 
 
@@ -233,6 +246,7 @@ Depending on the ``max_iterations`` parameter set, very high quality solutions c
         distance_matrix: np.ndarray,
         x0: Optional[List[int]] = None,
         max_iterations: Optional[int] = None,
+        rng: Optional[Random] = None,
         log_file: Optional[str] = None,
         verbose: bool = False,
     )
@@ -248,13 +262,18 @@ Depending on the ``max_iterations`` parameter set, very high quality solutions c
     x0
         Initial permutation. If not provided, it starts with a random path.
 
-    max_iterations
+    max_iterations {None}
         The maximum number of iterations for the algorithm. If not specified,
         it defaults to the number of nodes in the distance matrix.
 
-    log_file
+    rng {None}
+        A python random number generator which can be seeded to ensure 
+        consistent results given a call to this function. If not provided
+        the initial random number generator is used.
+
+    log_file {None}
         If not `None`, creates a log file with details about the whole
         execution.
 
-    verbose
+    verbose {False}
         If true, prints algorithm status every iteration.

python_tsp/heuristics/local_search.py

@@ -2,6 +2,7 @@
 
 from timeit import default_timer
 from typing import List, Optional, Tuple, TextIO
+from random import Random
 
 import numpy as np
 
@@ -20,6 +21,8 @@ def solve_tsp_local_search(
     x0: Optional[List[int]] = None,
     perturbation_scheme: str = "two_opt",
     max_processing_time: Optional[float] = None,
+    max_iterations: Optional[int] = None,
+    rng: Optional[Random] = None,
     log_file: Optional[str] = None,
     verbose: bool = False,
 ) -> Tuple[List, float]:
@@ -41,6 +44,10 @@ def solve_tsp_local_search(
         Maximum processing time in seconds. If not provided, the method stops
         only when a local minimum is obtained
 
+    rng
+        Random number generator to be passed to the pertubation scheme. If not
+        provided, the initial random generator is used.
+
     log_file
         If not `None`, creates a log file with details about the whole
         execution
@@ -79,7 +86,9 @@ def solve_tsp_local_search(
 
     while improvement and (not stop_early):
         improvement = False
-        for n_index, xn in enumerate(neighborhood_gen[perturbation_scheme](x)):
+        for n_index, xn in enumerate(
+            neighborhood_gen[perturbation_scheme](x, rng)
+        ):
             if default_timer() - tic > max_processing_time:
                 _print_message(TIME_LIMIT_MSG, verbose, log_file_handler)
                 stop_early = True

python_tsp/heuristics/perturbation_schemes.py

@@ -13,64 +13,78 @@
 .. [2] 2-opt: https://en.wikipedia.org/wiki/2-opt
 """
 
-from random import sample
-from typing import Callable, Dict, Generator, List
+from random import Random
+from typing import Callable, Dict, Generator, List, Union
 
+initial_random_generator = Random()
 
-def ps1_gen(x: List[int]) -> Generator[List[int], List[int], None]:
+
+def ps1_gen(
+    x: List[int], rng: Union[Random, None] = None
+) -> Generator[List[int], List[int], None]:
     """PS1 perturbation scheme: Swap two adjacent terms [1]
     This scheme has at most n - 1 swaps.
     """
+    rng = rng or initial_random_generator
 
     n = len(x)
     i_range = range(1, n - 1)
-    for i in sample(i_range, len(i_range)):
+    for i in rng.sample(i_range, len(i_range)):
         xn = x.copy()
         xn[i], xn[i + 1] = xn[i + 1], xn[i]
         yield xn
 
 
-def ps2_gen(x: List[int]) -> Generator[List[int], List[int], None]:
+def ps2_gen(
+    x: List[int], rng: Union[Random, None] = None
+) -> Generator[List[int], List[int], None]:
     """PS2 perturbation scheme: Swap any two elements [1]
     This scheme has n * (n - 1) / 2 swaps.
     """
+    rng = rng or initial_random_generator
 
     n = len(x)
     i_range = range(1, n - 1)
-    for i in sample(i_range, len(i_range)):
+    for i in rng.sample(i_range, len(i_range)):
         j_range = range(i + 1, n)
-        for j in sample(j_range, len(j_range)):
+        for j in rng.sample(j_range, len(j_range)):
             xn = x.copy()
             xn[i], xn[j] = xn[j], xn[i]
             yield xn
 
 
-def ps3_gen(x: List[int]) -> Generator[List[int], List[int], None]:
+def ps3_gen(
+    x: List[int], rng: Union[Random, None] = None
+) -> Generator[List[int], List[int], None]:
     """PS3 perturbation scheme: A single term is moved [1]
     This scheme has n * (n - 1) swaps.
     """
+    rng = rng or initial_random_generator
 
     n = len(x)
     i_range = range(1, n)
-    for i in sample(i_range, len(i_range)):
+    for i in rng.sample(i_range, len(i_range)):
         j_range = [j for j in range(1, n) if j != i]
-        for j in sample(j_range, len(j_range)):
+        for j in rng.sample(j_range, len(j_range)):
             xn = x.copy()
             node = xn.pop(i)
             xn.insert(j, node)
             yield xn
 
 
-def ps4_gen(x: List[int]) -> Generator[List[int], List[int], None]:
+def ps4_gen(
+    x: List[int], rng: Union[Random, None] = None
+) -> Generator[List[int], List[int], None]:
     """PS4 perturbation scheme: A subsequence is moved [1]"""
+    rng = rng or initial_random_generator
 
     n = len(x)
     i_range = range(1, n)
-    for i in sample(i_range, len(i_range)):
+    for i in rng.sample(i_range, len(i_range)):
         j_range = range(i + 1, n + 1)
-        for j in sample(j_range, len(j_range)):
+        for j in rng.sample(j_range, len(j_range)):
             k_range = [k for k in range(1, n + 1) if k not in range(i, j + 1)]
-            for k in sample(k_range, len(k_range)):
+            for k in rng.sample(k_range, len(k_range)):
                 xn = x.copy()
                 if k < i:
                     xn = xn[:k] + xn[i:j] + xn[k:i] + xn[j:]
@@ -79,29 +93,35 @@ def ps4_gen(x: List[int]) -> Generator[List[int], List[int], None]:
                 yield xn
 
 
-def ps5_gen(x: List[int]) -> Generator[List[int], List[int], None]:
+def ps5_gen(
+    x: List[int], rng: Union[Random, None] = None
+) -> Generator[List[int], List[int], None]:
     """PS5 perturbation scheme: A subsequence is reversed [1]"""
+    rng = rng or initial_random_generator
 
     n = len(x)
     i_range = range(1, n)
-    for i in sample(i_range, len(i_range)):
+    for i in rng.sample(i_range, len(i_range)):
         j_range = range(i + 2, n + 1)
-        for j in sample(j_range, len(j_range)):
+        for j in rng.sample(j_range, len(j_range)):
             xn = x.copy()
             xn = xn[:i] + list(reversed(xn[i:j])) + xn[j:]
             yield xn
 
 
-def ps6_gen(x: List[int]) -> Generator[List[int], List[int], None]:
+def ps6_gen(
+    x: List[int], rng: Union[Random, None] = None
+) -> Generator[List[int], List[int], None]:
     """PS6 perturbation scheme: A subsequence is reversed and moved [1]"""
+    rng = rng or initial_random_generator
 
     n = len(x)
     i_range = range(1, n)
-    for i in sample(i_range, len(i_range)):
+    for i in rng.sample(i_range, len(i_range)):
         j_range = range(i + 1, n + 1)
-        for j in sample(j_range, len(j_range)):
+        for j in rng.sample(j_range, len(j_range)):
             k_range = [k for k in range(1, n + 1) if k not in range(i, j + 1)]
-            for k in sample(k_range, len(k_range)):
+            for k in rng.sample(k_range, len(k_range)):
                 xn = x.copy()
                 if k < i:
                     xn = xn[:k] + list(reversed(xn[i:j])) + xn[k:i] + xn[j:]
@@ -110,21 +130,28 @@ def ps6_gen(x: List[int]) -> Generator[List[int], List[int], None]:
                 yield xn
 
 
-def two_opt_gen(x: List[int]) -> Generator[List[int], List[int], None]:
+def two_opt_gen(
+    x: List[int], rng: Union[Random, None] = None
+) -> Generator[List[int], List[int], None]:
     """2-opt perturbation scheme [2]"""
+    rng = rng or initial_random_generator
+
     n = len(x)
     i_range = range(2, n)
-    for i in sample(i_range, len(i_range)):
+    for i in rng.sample(i_range, len(i_range)):
         j_range = range(i + 1, n + 1)
-        for j in sample(j_range, len(j_range)):
+        for j in rng.sample(j_range, len(j_range)):
             xn = x.copy()
             xn = xn[: i - 1] + list(reversed(xn[i - 1 : j])) + xn[j:]
             yield xn
 
 
 # Mapping with all possible neighborhood generators in this module
 neighborhood_gen: Dict[
-    str, Callable[[List[int]], Generator[List[int], List[int], None]]
+    str,
+    Callable[
+        [List[int], Union[Random, None]], Generator[List[int], List[int], None]
+    ],
 ] = {
     "ps1": ps1_gen,
     "ps2": ps2_gen,

python_tsp/heuristics/record_to_record.py

@@ -1,4 +1,4 @@
-from random import randint
+from random import Random, randint
 from typing import List, Optional, TextIO
 
 import numpy as np
@@ -21,6 +21,7 @@ def solve_tsp_record_to_record(
     distance_matrix: np.ndarray,
     x0: Optional[List[int]] = None,
     max_iterations: Optional[int] = None,
+    rng: Optional[Random] = None,
     log_file: Optional[str] = None,
     verbose: bool = False,
 ):
@@ -41,6 +42,10 @@ def solve_tsp_record_to_record(
         The maximum number of iterations for the algorithm. If not specified,
         it defaults to the number of nodes in the distance matrix.
 
+    rng
+        Random number generator to be passed to the pertubation scheme. If not
+        provided, the initial random generator is used.
+
     log_file
         If not `None`, creates a log file with details about the whole
         execution.
@@ -62,15 +67,17 @@ def solve_tsp_record_to_record(
     max_iterations = max_iterations or n
     x, fx = setup_initial_solution(distance_matrix=distance_matrix, x0=x0)
 
+    new_int = rng.randint if rng else randint
+
     log_file_handler = (
         open(log_file, "w", encoding="utf-8") if log_file else None
     )
 
     for iteration in range(1, max_iterations + 1):
         xn = x[:]
         for _ in range(2):
-            u = randint(1, n - 1)
-            v = randint(1, n - 1)
+            u = new_int(1, n - 1)
+            v = new_int(1, n - 1)
             xn[u], xn[v] = xn[v], xn[u]
 
         xn, fn = solve_tsp_lin_kernighan(