Operator dependencies (#70)

* Add tests for operator dependencies

* Implement operator coupling

* Improve docstrings and variable names

alns/ALNS.py

@@ -1,7 +1,9 @@
+from collections import defaultdict
 import logging
 import time
-from typing import Callable, Dict, List, Optional, Tuple
+from typing import Callable, Dict, Iterable, List, Optional, Tuple
 
+import numpy as np
 import numpy.random as rnd
 
 from alns.Result import Result
@@ -55,6 +57,8 @@ def __init__(self, rnd_state: rnd.RandomState = rnd.RandomState()):
 
         self._rnd_state = rnd_state
 
+        self._only_after: Dict[_OperatorType, set] = defaultdict(set)
+
     @property
     def destroy_operators(self) -> List[Tuple[str, _OperatorType]]:
         """
@@ -96,7 +100,13 @@ def add_destroy_operator(self, op: _OperatorType, name: str = None):
         logger.debug(f"Adding destroy operator {op.__name__}.")
         self._destroy_operators[op.__name__ if name is None else name] = op
 
-    def add_repair_operator(self, op: _OperatorType, name: str = None):
+    def add_repair_operator(
+        self,
+        op: _OperatorType,
+        name: str = None,
+        *,
+        only_after: Optional[Iterable[_OperatorType]] = None,
+    ):
         """
         Adds a repair operator to the heuristic instance.
 
@@ -109,10 +119,49 @@ def add_repair_operator(self, op: _OperatorType, name: str = None):
         name
             Optional name argument, naming the operator. When not passed, the
             function name is used instead.
+        only_after
+            Optional keyword-only argument indicating which destroy operators
+            work with the passed-in repair operator. If passed, this argument
+            should be an iterable (e.g. a list) of destroy operators. If not
+            passed, the default is to assume that all destroy operators work
+            with the new repair operator.
         """
         logger.debug(f"Adding repair operator {op.__name__}.")
         self._repair_operators[name if name else op.__name__] = op
 
+        if only_after is not None:
+            self._only_after[op].update(only_after)
+
+    def _compute_op_coupling(self) -> np.ndarray:
+        """
+        Internal helper to compute a matrix that describes the
+        coupling between destroy and repair operators. The matrix has size
+        |d_ops|-by-|r_ops| and entry (i, j) is 1 if destroy operator i can
+        be used in conjunction with repair operator j and 0 otherwise.
+
+        If the only_after keyword-only argument was not used when adding
+        the repair operators, then all entries of the matrix are 1.
+        """
+        op_coupling = np.ones(
+            (len(self.destroy_operators), len(self.repair_operators))
+        )
+
+        for r_idx, (_, r_op) in enumerate(self.repair_operators):
+            coupled_d_ops = self._only_after[r_op]
+
+            for d_idx, (_, d_op) in enumerate(self.destroy_operators):
+                if coupled_d_ops and d_op not in coupled_d_ops:
+                    op_coupling[d_idx, r_idx] = 0
+
+        # Destroy operators must be coupled with at least one repair operator
+        d_idcs = np.flatnonzero(np.count_nonzero(op_coupling, axis=1) == 0)
+
+        if d_idcs.size != 0:
+            d_name, _ = self.destroy_operators[d_idcs[0]]
+            raise ValueError(f"{d_name} has no coupled repair operators.")
+
+        return op_coupling
+
     def iterate(
         self,
         initial_solution: State,
@@ -170,6 +219,7 @@ class of vehicle routing problems with backhauls. *European Journal of
 
         curr = best = initial_solution
         init_obj = initial_solution.objective()
+        op_coupling = self._compute_op_coupling()
 
         logger.debug(f"Initial solution has objective {init_obj:.2f}.")
 
@@ -178,7 +228,9 @@ class of vehicle routing problems with backhauls. *European Journal of
         stats.collect_runtime(time.perf_counter())
 
         while not stop(self._rnd_state, best, curr):
-            d_idx, r_idx = weight_scheme.select_operators(self._rnd_state)
+            d_idx, r_idx = weight_scheme.select_operators(
+                self._rnd_state, op_coupling
+            )
 
             d_name, d_operator = self.destroy_operators[d_idx]
             r_name, r_operator = self.repair_operators[r_idx]

alns/tests/test_alns.py

@@ -1,3 +1,4 @@
+import numpy as np
 import numpy.random as rnd
 from numpy.testing import (
     assert_,
@@ -48,6 +49,34 @@ def objective(self):
         return self._value
 
 
+def get_repair_operators(n):
+    """
+    Get a list of n dummy repair operators.
+    """
+    repair_operators = []
+
+    for idx in range(n):
+        op = lambda: None
+        op.__name__ = f"Repair operator {idx}"
+        repair_operators.append(op)
+
+    return repair_operators
+
+
+def get_destroy_operators(n):
+    """
+    Get a list of n dummy destroy operators.
+    """
+    destroy_operators = []
+
+    for idx in range(n):
+        op = lambda: None
+        op.__name__ = f"Destroy operator {idx}"
+        destroy_operators.append(op)
+
+    return destroy_operators
+
+
 # CALLBACKS --------------------------------------------------------------------
 
 
@@ -134,6 +163,68 @@ def repair_operator():  # placeholder
     assert_(operator is repair_operator)
 
 
+def test_compute_op_coupling():
+    """
+    Tests if the compute_op_coupling method correctly computes the matrix
+    that describes the dependencies between repair and destroy operators.
+    """
+    alns = ALNS()
+
+    d_operators = get_destroy_operators(2)
+    r_operators = get_repair_operators(2)
+
+    for d_op in d_operators:
+        alns.add_destroy_operator(d_op)
+
+    for r_op in r_operators:
+        alns.add_repair_operator(r_op)
+
+    op_coupling = alns._compute_op_coupling()
+
+    assert_almost_equal(op_coupling, np.ones((2, 2)))
+
+
+def test_compute_op_coupling_only_after():
+    """
+    Tests if the compute_op_coupling method correctly computes the matrix
+    when the only_after paramter is specified for certain repair operators.
+    """
+    alns = ALNS()
+
+    d_operators = get_destroy_operators(2)
+    r_operators = get_repair_operators(2)
+
+    for d_op in d_operators:
+        alns.add_destroy_operator(d_op)
+
+    for idx, r_op in enumerate(r_operators):
+        alns.add_repair_operator(r_op, only_after=[d_operators[idx]])
+
+    op_coupling = alns._compute_op_coupling()
+
+    assert_almost_equal(op_coupling, np.eye(2))
+
+
+def test_raise_uncoupled_destroy_op():
+    """
+    Tests if having a destroy operator that is not coupled to any of the
+    repair operators raises an an error when computing the operator coupling.
+    """
+    alns = ALNS()
+
+    d_operators = get_destroy_operators(2)
+    r_operators = get_repair_operators(2)
+
+    for d_op in d_operators:
+        alns.add_destroy_operator(d_op)
+
+    for r_op in r_operators:
+        alns.add_repair_operator(r_op, only_after=[d_operators[0]])
+
+    with assert_raises(ValueError):
+        alns._compute_op_coupling()
+
+
 # PARAMETERS -------------------------------------------------------------------
 
 

alns/weights/SegmentedWeights.py

@@ -54,7 +54,7 @@ def __init__(
 
         self._reset_segment_weights()
 
-    def select_operators(self, rnd_state):
+    def select_operators(self, rnd_state, op_coupling):
         self._iter += 1
 
         if self._iter % self._seg_length == 0:
@@ -68,7 +68,7 @@ def select_operators(self, rnd_state):
 
             self._reset_segment_weights()
 
-        return super().select_operators(rnd_state)
+        return super().select_operators(rnd_state, op_coupling)
 
     def update_weights(self, d_idx, r_idx, s_idx):
         self._d_seg_weights[d_idx] += self._scores[s_idx]

alns/weights/WeightScheme.py

@@ -37,7 +37,9 @@ def destroy_weights(self) -> np.ndarray:
     def repair_weights(self) -> np.ndarray:
         return self._r_weights
 
-    def select_operators(self, rnd_state: RandomState) -> Tuple[int, int]:
+    def select_operators(
+        self, rnd_state: RandomState, op_coupling: np.ndarray
+    ) -> Tuple[int, int]:
         """
         Selects a destroy and repair operator pair to apply in this iteration.
         The default implementation uses a roulette wheel mechanism, where each
@@ -48,6 +50,11 @@ def select_operators(self, rnd_state: RandomState) -> Tuple[int, int]:
         rnd_state
             Random state object, to be used for number generation.
 
+        op_coupling
+            Matrix that indicates coupling between destroy and repair
+            operators. Entry (i, j) is 1 if destroy operator i can be used in
+            conjunction with repair operator j and 0 otherwise.
+
         Returns
         -------
         A tuple of (d_idx, r_idx), which are indices into the destroy and repair
@@ -58,7 +65,11 @@ def select(op_weights):
             probs = op_weights / np.sum(op_weights)
             return rnd_state.choice(range(len(op_weights)), p=probs)
 
-        return select(self._d_weights), select(self._r_weights)
+        d_idx = select(self._d_weights)
+        coupled_r_idcs = np.flatnonzero(op_coupling[d_idx])
+        r_idx = coupled_r_idcs[select(self._r_weights[coupled_r_idcs])]
+
+        return d_idx, r_idx
 
     @abstractmethod
     def update_weights(self, d_idx: int, r_idx: int, s_idx: int):

alns/weights/tests/test_segmented_weights.py

@@ -34,10 +34,11 @@ def test_update_weights(
 ):
     rnd_state = np.random.RandomState(1)
     weights = SegmentedWeights(scores, 1, 1, seg_decay, 1)
+    op_coupling = np.ones((1, 1))
 
     # TODO other weights?
     weights.update_weights(0, 0, 1)
-    weights.select_operators(rnd_state)
+    weights.select_operators(rnd_state, op_coupling)
 
     assert_almost_equal(weights.destroy_weights[0], expected[0])
     assert_almost_equal(weights.repair_weights[0], expected[1])

alns/weights/tests/test_simple_weights.py

@@ -1,7 +1,8 @@
 from typing import List
 
 import numpy as np
-from numpy.testing import assert_raises, assert_almost_equal
+import numpy.random as rnd
+from numpy.testing import assert_, assert_raises, assert_almost_equal
 from pytest import mark
 
 from alns.weights import SimpleWeights
@@ -38,4 +39,25 @@ def test_update_weights(
     assert_almost_equal(weights.repair_weights[0], expected[1])
 
 
-# TODO test select weights
+@mark.parametrize(
+    "op_coupling",
+    [
+        np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]]),
+        np.array([[1, 0, 1], [0, 1, 1], [1, 1, 0]]),
+        np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+        np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]]),
+        np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]]),
+        np.array([[1, 0, 0], [1, 0, 0], [1, 0, 0]]),  # Not allowed by ALNS
+    ],
+)
+def test_select_operators(op_coupling):
+    """
+    Test if the indices of the selected operators correspond to the
+    ones that are given by the operator coupling.
+    """
+    rnd_state = rnd.RandomState()
+    n_destroy, n_repair = op_coupling.shape
+    weights = SimpleWeights([0, 0, 0, 0], n_destroy, n_repair, 0)
+    d_idx, r_idx = weights.select_operators(rnd_state, op_coupling)
+
+    assert_((d_idx, r_idx) in np.argwhere(op_coupling == 1))