Implement TSP generator and LazyConstraintsComponent

miplearn/__init__.py

@@ -10,6 +10,7 @@ from .extractors import (SolutionExtractor,
                         )
 from .components.component import Component
 from .components.objective import ObjectiveValueComponent
+from .components.lazy import LazyConstraintsComponent
 from .components.primal import (PrimalSolutionComponent,
                                 AdaptivePredictor,
                                )

miplearn/components/branching.py

@@ -118,29 +118,3 @@ class BranchPriorityComponent(Component):
         instance_features = instance.get_instance_features()
         var_features = instance.get_variable_features(var, index)
         return np.hstack([instance_features, var_features])
-            
-    def merge(self, other_components):
-        keys = set(self.x_train.keys())
-        for comp in other_components:
-            self.pending_instances += comp.pending_instances
-            keys = keys.union(set(comp.x_train.keys()))
-            
-        # Merge x_train and y_train
-        for key in keys:
-            x_train_submatrices = [comp.x_train[key]
-                                   for comp in other_components
-                                   if key in comp.x_train.keys()]
-            y_train_submatrices = [comp.y_train[key]
-                                   for comp in other_components
-                                   if key in comp.y_train.keys()]
-            if key in self.x_train.keys():
-                x_train_submatrices += [self.x_train[key]]
-                y_train_submatrices += [self.y_train[key]]
-            self.x_train[key] = np.vstack(x_train_submatrices)
-            self.y_train[key] = np.vstack(y_train_submatrices)
-            
-        # Merge trained ML predictors
-        for comp in other_components:
-            for key in comp.predictors.keys():
-                if key not in self.predictors.keys():
-                    self.predictors[key] = comp.predictors[key]

miplearn/components/component.py

@@ -18,10 +18,6 @@ class Component(ABC):
     def after_solve(self, solver, instance, model):
         pass
     
-    @abstractmethod
-    def merge(self, other):
-        pass
-    
     @abstractmethod
     def fit(self, training_instances):
         pass

miplearn/components/lazy.py

@@ -0,0 +1,48 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+from .component import Component
+from ..extractors import *
+
+from abc import ABC, abstractmethod
+from copy import deepcopy
+import numpy as np
+from sklearn.pipeline import make_pipeline
+from sklearn.linear_model import LogisticRegression
+from sklearn.preprocessing import StandardScaler
+from sklearn.model_selection import cross_val_score
+from sklearn.metrics import roc_curve
+from sklearn.neighbors import KNeighborsClassifier
+from tqdm.auto import tqdm
+import pyomo.environ as pe
+import logging
+logger = logging.getLogger(__name__)
+
+
+class LazyConstraintsComponent(Component):
+    """
+    A component that predicts which lazy constraints to enforce.
+    """
+    
+    def __init__(self):
+        self.violations = set()
+    
+    def before_solve(self, solver, instance, model):
+        logger.info("Enforcing %d lazy constraints" % len(self.violations))
+        for v in self.violations:
+            cut = instance.build_lazy_constraint(model, v)
+            solver.internal_solver.add_constraint(cut)
+        
+    def after_solve(self, solver, instance, model):
+        pass
+                
+    def fit(self, training_instances):
+        for instance in training_instances:
+            if not hasattr(instance, "found_violations"):
+                continue
+            for v in instance.found_violations:
+                self.violations.add(v)
+                
+    def predict(self, instance, model=None):
+        return self.violations

miplearn/components/objective.py

@@ -30,9 +30,6 @@ class ObjectiveValueComponent(Component):
     def after_solve(self, solver, instance, model):
         pass
     
-    def merge(self, other):
-        pass
-    
     def fit(self, training_instances):
         features = InstanceFeaturesExtractor().extract(training_instances)
         ub = ObjectiveValueExtractor(kind="upper bound").extract(training_instances)

miplearn/components/primal.py

@@ -200,6 +200,3 @@ class PrimalSolutionComponent(Component):
                     if ws[i, 1] >= self.thresholds[category, label]:
                         solution[var][index] = label
         return solution
-
-    def merge(self, other_components):
-        pass

miplearn/problems/tests/test_tsp.py

@@ -0,0 +1,68 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+from miplearn import LearningSolver
+from miplearn.problems.tsp import TravelingSalesmanGenerator, TravelingSalesmanInstance
+import numpy as np
+from numpy.linalg import norm
+from scipy.spatial.distance import pdist, squareform
+from scipy.stats import uniform, randint
+
+
+def test_generator():
+    instances = TravelingSalesmanGenerator(x=uniform(loc=0.0, scale=1000.0),
+                                           y=uniform(loc=0.0, scale=1000.0),
+                                           n=randint(low=100, high=101),
+                                           gamma=uniform(loc=0.95, scale=0.1),
+                                           fix_cities=True).generate(100)
+    assert len(instances) == 100
+    assert instances[0].n_cities == 100
+    assert norm(instances[0].distances - instances[0].distances.T) < 1e-6
+    d = [instance.distances[0,1] for instance in instances]
+    assert np.std(d) > 0
+
+    
+def test_instance():
+    n_cities = 4
+    distances = np.array([
+        [0., 1., 2., 1.],
+        [1., 0., 1., 2.],
+        [2., 1., 0., 1.],
+        [1., 2., 1., 0.],
+    ])
+    instance = TravelingSalesmanInstance(n_cities, distances)
+    solver = LearningSolver()
+    solver.solve(instance)
+    x = instance.solution["x"]
+    assert x[0,1] == 1.0
+    assert x[0,2] == 0.0
+    assert x[0,3] == 1.0
+    assert x[1,2] == 1.0
+    assert x[1,3] == 0.0
+    assert x[2,3] == 1.0
+    assert instance.lower_bound == 4.0
+    assert instance.upper_bound == 4.0
+
+    
+def test_subtour():
+    n_cities = 6
+    cities = np.array([
+        [0., 0.],
+        [1., 0.],
+        [2., 0.],
+        [3., 0.],
+        [0., 1.],
+        [3., 1.],
+    ])
+    distances = squareform(pdist(cities))
+    instance = TravelingSalesmanInstance(n_cities, distances)
+    solver = LearningSolver()
+    solver.solve(instance)
+    x = instance.solution["x"]
+    assert x[0,1] == 1.0
+    assert x[0,4] == 1.0
+    assert x[1,2] == 1.0
+    assert x[2,3] == 1.0
+    assert x[3,5] == 1.0
+    assert x[4,5] == 1.0

miplearn/problems/tsp.py

@@ -5,71 +5,109 @@
 import numpy as np
 import pyomo.environ as pe
 from miplearn import Instance
+from scipy.stats import uniform, randint
+from scipy.spatial.distance import pdist, squareform
+from scipy.stats.distributions import rv_frozen
+import networkx as nx
 import random
 
 
-class TravelingSalesmanChallengeA:
-    """Fixed set of cities, small perturbation to travel speed."""
-    def __init__():
-        self.generator = TravelingSalesmanGenerator(speed=uniform(loc=0.9, scale=0.2),
-                                                    x=uniform(loc=0.0, loc=1000.0),
-                                                    y=uniform(loc=0.0, loc=1000.0),
-                                                    pn=0.0,
-                                                    n=randint(low=100, high=100),
-                                                    fix_cities=True)
+class ChallengeA:
+    def __init__(self,
+                 seed=42,
+                 n_training_instances=500,
+                 n_test_instances=50,
+                ):
         
-    def get_training_instances():
-        return self.generator.generate(500)
-    
-    def get_test_instances():
-        return self.generator.generate(100)
+        np.random.seed(seed)
+        self.generator = TravelingSalesmanGenerator(x=uniform(loc=0.0, scale=1000.0),
+                                                    y=uniform(loc=0.0, scale=1000.0),
+                                                    n=randint(low=350, high=351),
+                                                    gamma=uniform(loc=0.95, scale=0.1),
+                                                    fix_cities=True,
+                                                   )
+
+        np.random.seed(seed + 1)
+        self.training_instances = self.generator.generate(n_training_instances)
+        
+        np.random.seed(seed + 2)
+        self.test_instances = self.generator.generate(n_test_instances)        
     
         
 class TravelingSalesmanGenerator:
-    """Random generator for the Traveling Salesman Problem.
+    """Random generator for the Traveling Salesman Problem."""
 
-    The generator starts by randomly selecing n points with coordinates (x_i, y_i), where n, x_i 
-    and y_i are random variables. The time required to travel from a pair of cities is calculated
-    by: (i) computing the euclidean distance between the cities, (ii) sampling a random variable
-    speed_i, (iii) dividing the two numbers.
+    def __init__(self,
+                 x=uniform(loc=0.0, scale=1000.0),
+                 y=uniform(loc=0.0, scale=1000.0),
+                 n=randint(low=100, high=101),
+                 gamma=uniform(loc=1.0, scale=0.0),
+                 fix_cities=True,
+                 round=True,
+                ):
+        """Initializes the problem generator.
+        
+        Initially, the generator creates n cities (x_1,y_1),...,(x_n,y_n) where n, x_i and y_i are
+        sampled independently from the provided probability distributions `n`, `x` and `y`. For each
+        (unordered) pair of cities (i,j), the distance d[i,j] between them is set to:
 
-    If fix_cities is True, the cities and travel times will be calculated only once, during the
-    constructor. Each time an instance is generated, however, each city will have probability pv
-    of being removed from the list. If fix_cities if False, then the cities and travel times will
-    be resampled each time an instance is generated. The probability pn is not used in this case.
+            d[i,j] = gamma[i,j] \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2}
 
-    All random variables are independent.
-    """
+        where gamma is sampled from the provided probability distribution `gamma`.
 
-    def __init__(self,
-                 speed=uniform(loc=0.75, scale=0.5),
-                 x=uniform(loc=0.0, loc=1000.0),
-                 y=uniform(loc=0.0, loc=1000.0),
-                 pn=0.0,
-                 n=randint(low=100, high=100),
-                 fix_cities=True):
-        """Initializes the problem generator.
+        If fix_cities=True, the list of cities is kept the same for all generated instances. The
+        gamma values, and therefore also the distances, are still different.
+        
+        By default, all distances d[i,j] are rounded to the nearest integer.  If `round=False`
+        is provided, this rounding will be disabled.
 
         Arguments
         ---------
-        speed: rv_continuous
-            Probability distribution for travel speed.
         x: rv_continuous
             Probability distribution for the x-coordinate of each city.
         y: rv_continuous
             Probability distribution for the y-coordinate of each city.
-        pn: float
-            Probability of a city being removed from the list. Only used if fix_cities=True.
         n: rv_discrete
             Probability distribution for the number of cities.
         fix_cities: bool
-            If true, cities will be resampled for every generated instance. Otherwise, list of
+            If False, cities will be resampled for every generated instance. Otherwise, list of
             cities will be computed once, during the constructor.
+        round: bool
+            If True, distances are rounded to the nearest integer.
         """
-        pass
+        assert isinstance(x, rv_frozen), "x should be a SciPy probability distribution"
+        assert isinstance(y, rv_frozen), "y should be a SciPy probability distribution"
+        assert isinstance(n, rv_frozen), "n should be a SciPy probability distribution"
+        assert isinstance(gamma, rv_frozen), "gamma should be a SciPy probability distribution"
+        self.x = x
+        self.y = y
+        self.n = n
+        self.gamma = gamma
+        self.round = round
+        
+        if fix_cities:
+            self.fixed_n, self.fixed_cities = self._generate_cities()
+        else:
+            self.fixed_n = None
+            self.fixed_cities = None
     
     def generate(self, n_samples):
-        pass
+        def _sample():
+            if self.fixed_cities is not None:
+                n, cities = self.fixed_n, self.fixed_cities
+            else:
+                n, cities = self._generate_cities()
+            distances = squareform(pdist(cities)) * self.gamma.rvs(size=(n, n))
+            distances = np.tril(distances) + np.triu(distances.T, 1)
+            if self.round:
+                distances = distances.round()
+            return TravelingSalesmanInstance(n, distances)
+        return [_sample() for _ in range(n_samples)]
+            
+    def _generate_cities(self):
+        n = self.n.rvs()
+        cities = np.array([(self.x.rvs(), self.y.rvs()) for _ in range(n)])
+        return n, cities
     
 
 class TravelingSalesmanInstance(Instance):
@@ -82,9 +120,49 @@ class TravelingSalesmanInstance(Instance):
     """
     
     def __init__(self, n_cities, distances):
-        assert isinstance(distances, np.array)
+        assert isinstance(distances, np.ndarray)
         assert distances.shape == (n_cities, n_cities)
         self.n_cities = n_cities
         self.distances = distances
-        pass
     
+    def to_model(self):
+        self.model = model = pe.ConcreteModel()
+        self.edges = edges = [(i,j)
+                              for i in range(self.n_cities)
+                              for j in range(i+1, self.n_cities)]
+        model.x = pe.Var(edges, domain=pe.Binary)
+        model.obj = pe.Objective(rule=lambda m : sum(m.x[i,j] * self.distances[i,j]
+                                                     for (i,j) in edges),
+                                 sense=pe.minimize)
+        model.eq_degree = pe.ConstraintList()
+        model.eq_subtour = pe.ConstraintList()
+        for i in range(self.n_cities):
+            model.eq_degree.add(sum(model.x[min(i,j), max(i,j)]
+                                    for j in range(self.n_cities) if i != j) == 2)
+        return model
+    
+    def get_instance_features(self):
+        return np.array([1])
+    
+    def get_variable_features(self, var, index):
+        return np.array([1])
+    
+    def get_variable_category(self, var, index):
+        return index
+    
+    def find_violations(self, model):
+        selected_edges = [e for e in self.edges if model.x[e].value > 0.5]
+        graph = nx.Graph()
+        graph.add_edges_from(selected_edges)
+        components = [frozenset(c) for c in list(nx.connected_components(graph))]
+        violations = []
+        for c in components:
+            if len(c) < self.n_cities:
+                violations += [c]
+        return violations
+    
+    def build_lazy_constraint(self, model, component):
+        cut_edges = [e for e in self.edges
+                     if (e[0] in component and e[1] not in component) or
+                        (e[0] not in component and e[1] in component)]
+        return model.eq_subtour.add(sum(model.x[e] for e in cut_edges) >= 2)

miplearn/solvers.py

@@ -2,7 +2,7 @@
 #  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 #  Released under the modified BSD license. See COPYING.md for more details.
 
-from . import ObjectiveValueComponent, PrimalSolutionComponent
+from . import ObjectiveValueComponent, PrimalSolutionComponent, LazyConstraintsComponent
 import pyomo.environ as pe
 from pyomo.core import Var
 from copy import deepcopy
@@ -88,6 +88,9 @@ class InternalSolver:
                 self.solver.update_var(var[index])        
         logger.info("Fixing values for %d variables (out of %d)" %
                     (count_fixed, count_total))
+        
+    def add_constraint(self, cut):
+        self.solver.add_constraint(cut)
 
     
 class GurobiSolver(InternalSolver):
@@ -198,6 +201,7 @@ class LearningSolver:
             self.components = {
                 "ObjectiveValue": ObjectiveValueComponent(),
                 "PrimalSolution": PrimalSolutionComponent(),
+                "LazyConstraints": LazyConstraintsComponent(),
             }
             
         assert self.mode in ["exact", "heuristic"]
@@ -231,27 +235,44 @@ class LearningSolver:
         self.internal_solver = self._create_internal_solver()
         self.internal_solver.set_model(model)
 
-        # Solve LP relaxation
+        logger.debug("Solving LP relaxation...")
         results = self.internal_solver.solve_lp(tee=tee)
         instance.lp_solution = self.internal_solver.get_solution()
         instance.lp_value = results["Optimal value"]
         
-        # Invoke before_solve callbacks
+        logger.debug("Running before_solve callbacks...")
         for component in self.components.values():
             component.before_solve(self, instance, model)
         
         if relaxation_only:
             return results
         
-        # Solver original MIP
-        results = self.internal_solver.solve(tee=tee)
+        total_wallclock_time = 0
+        instance.found_violations = []
+        while True:
+            logger.debug("Solving MIP...")
+            results = self.internal_solver.solve(tee=tee)
+            logger.debug("    %.2f s" % results["Wallclock time"])
+            total_wallclock_time += results["Wallclock time"]
+            if not hasattr(instance, "find_violations"):
+                break
+            logger.debug("Finding violated constraints...")    
+            violations = instance.find_violations(model)
+            if len(violations) == 0:
+                break
+            instance.found_violations += violations
+            logger.debug("    %d violations found" % len(violations))
+            for v in violations:
+                cut = instance.build_lazy_constraint(model, v)
+                self.internal_solver.add_constraint(cut)
+        results["Wallclock time"] = total_wallclock_time
         
         # Read MIP solution and bounds
         instance.lower_bound = results["Lower bound"]
         instance.upper_bound = results["Upper bound"]
         instance.solution = self.internal_solver.get_solution()
         
-        # Invoke after_solve callbacks
+        logger.debug("Calling after_solve callbacks...")    
         for component in self.components.values():
             component.after_solve(self, instance, model)
             
@@ -282,6 +303,7 @@ class LearningSolver:
                 "LP value": instance.lp_value,
                 "Upper bound": instance.upper_bound,
                 "Lower bound": instance.lower_bound,
+                "Violations": instance.found_violations,
             }
 
         p_map_results = p_map(_process, instances, num_cpus=n_jobs, desc=label)
@@ -294,12 +316,7 @@ class LearningSolver:
             instances[idx].lp_value = r["LP value"]
             instances[idx].lower_bound = r["Lower bound"]
             instances[idx].upper_bound = r["Upper bound"]
-        
-        for (name, component) in self.components.items():
-            subcomponents = [subsolver.components[name]
-                             for subsolver in subsolvers
-                             if name in subsolver.components.keys()]
-            self.components[name].merge(subcomponents)
+            instances[idx].found_violations = r["Violations"]
         
         return results