diff --git a/.gitignore b/.gitignore index 330ca27..fd0123c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ +TODO.md +*.bin *$py.class *.cover *.egg diff --git a/miplearn/__init__.py b/miplearn/__init__.py index ab58ecc..48a470f 100644 --- a/miplearn/__init__.py +++ b/miplearn/__init__.py @@ -3,4 +3,4 @@ # Written by Alinson S. Xavier from .instance import Instance -from .solvers import LearningSolver \ No newline at end of file +from .solvers import LearningSolver diff --git a/miplearn/instance.py b/miplearn/instance.py index 49fa5d4..40c454c 100644 --- a/miplearn/instance.py +++ b/miplearn/instance.py @@ -4,23 +4,24 @@ from abc import ABC, abstractmethod + class Instance(ABC): """ Abstract class holding all the data necessary to generate a concrete model of the problem. In the knapsack problem, for example, this class could hold the number of items, their weights and costs, as well as the size of the knapsack. Objects implementing this class are able to - convert themselves into a concrete optimization model, which can be optimized by solver, or + convert themselves into a concrete optimization model, which can be optimized by a solver, or into arrays of features, which can be provided as inputs to machine learning models. """ - + @abstractmethod def to_model(self): """ Returns a concrete Pyomo model corresponding to this instance. """ pass - + @abstractmethod def get_instance_features(self): """ @@ -65,4 +66,4 @@ class Instance(ABC): pass def get_variable_category(self, var, index): - return "default" \ No newline at end of file + return "default" diff --git a/miplearn/problems/knapsack.py b/miplearn/problems/knapsack.py index 15271d3..eef184e 100644 --- a/miplearn/problems/knapsack.py +++ b/miplearn/problems/knapsack.py @@ -6,46 +6,48 @@ import miplearn import numpy as np import pyomo.environ as pe + class KnapsackInstance(miplearn.Instance): def __init__(self, weights, prices, capacity): self.weights = weights self.prices = prices self.capacity = capacity - + def to_model(self): - model = m = pe.ConcreteModel() + model = pe.ConcreteModel() items = range(len(self.weights)) - m.x = pe.Var(items, domain=pe.Binary) - m.OBJ = pe.Objective(rule=lambda m : sum(m.x[v] * self.prices[v] for v in items), - sense=pe.maximize) - m.eq_capacity = pe.Constraint(rule = lambda m : - sum(m.x[v] * self.weights[v] - for v in items) <= self.capacity) - return m - + model.x = pe.Var(items, domain=pe.Binary) + model.OBJ = pe.Objective(rule=lambda m: sum(m.x[v] * self.prices[v] for v in items), + sense=pe.maximize) + model.eq_capacity = pe.Constraint(rule=lambda m: sum(m.x[v] * self.weights[v] + for v in items) <= self.capacity) + return model + def get_instance_features(self): return np.array([ self.capacity, np.average(self.weights), ]) - + def get_variable_features(self, var, index): return np.array([ self.weights[index], self.prices[index], ]) - + + class KnapsackInstance2(KnapsackInstance): """ Alternative implementation of the Knapsack Problem, which assigns a different category for each decision variable, and therefore trains one machine learning model per variable. """ + def get_instance_features(self): return np.hstack([self.weights, self.prices]) - + def get_variable_features(self, var, index): return np.array([ ]) - + def get_variable_category(self, var, index): - return index \ No newline at end of file + return index diff --git a/miplearn/problems/stab.py b/miplearn/problems/stab.py deleted file mode 100644 index aa0e4b4..0000000 --- a/miplearn/problems/stab.py +++ /dev/null @@ -1,60 +0,0 @@ -# MIPLearn: A Machine-Learning Framework for Mixed-Integer Optimization -# Copyright (C) 2019-2020 Argonne National Laboratory. All rights reserved. -# Written by Alinson S. Xavier - -import numpy as np -import pyomo.environ as pe -import networkx as nx -from miplearn import Instance -import random - -class MaxStableSetGenerator: - def __init__(self, sizes=[50], densities=[0.1]): - self.sizes = sizes - self.densities = densities - - def generate(self): - size = random.choice(self.sizes) - density = random.choice(self.densities) - self.graph = nx.generators.random_graphs.binomial_graph(size, density) - weights = np.ones(self.graph.number_of_nodes()) - return MaxStableSetInstance(self.graph, weights) - - -class MaxStableSetInstance(Instance): - def __init__(self, graph, weights): - self.graph = graph - self.weights = weights - - def to_model(self): - nodes = list(self.graph.nodes) - edges = list(self.graph.edges) - model = m = pe.ConcreteModel() - m.x = pe.Var(nodes, domain=pe.Binary) - m.OBJ = pe.Objective(rule=lambda m : sum(m.x[v] * self.weights[v] for v in nodes), - sense=pe.maximize) - m.edge_eqs = pe.ConstraintList() - for edge in edges: - m.edge_eqs.add(m.x[edge[0]] + m.x[edge[1]] <= 1) - return m - - def get_instance_features(self): - return np.array([ - self.graph.number_of_nodes(), - self.graph.number_of_edges(), - ]) - - def get_variable_features(self, var, index): - first_neighbors = list(self.graph.neighbors(index)) - second_neighbors = [list(self.graph.neighbors(u)) for u in first_neighbors] - degree = len(first_neighbors) - neighbor_degrees = sorted([len(nn) for nn in second_neighbors]) - neighbor_degrees = neighbor_degrees + [100.] * 10 - return np.array([ - degree, - neighbor_degrees[0] - degree, - neighbor_degrees[1] - degree, - neighbor_degrees[2] - degree, - ]) - - diff --git a/miplearn/solvers.py b/miplearn/solvers.py index c783216..9e2bcb5 100644 --- a/miplearn/solvers.py +++ b/miplearn/solvers.py @@ -2,77 +2,81 @@ # Copyright (C) 2019-2020 Argonne National Laboratory. All rights reserved. # Written by Alinson S. Xavier -from .warmstart import * +# from .warmstart import WarmStartPredictor +from .transformers import PerVariableTransformer +from .warmstart import WarmStartPredictor import pyomo.environ as pe import numpy as np -from math import isfinite + class LearningSolver: """ Mixed-Integer Linear Programming (MIP) solver that extracts information from previous runs, using Machine Learning methods, to accelerate the solution of new (yet unseen) instances. """ - + def __init__(self, - threads = 4, - ws_predictor = None): - self.parent_solver = pe.SolverFactory('cplex_persistent') + threads=4, + parent_solver=pe.SolverFactory('cbc')): + self.parent_solver = parent_solver self.parent_solver.options["threads"] = threads - self.train_x = None - self.train_y = None - self.ws_predictor = ws_predictor - - def solve(self, - instance, - tee=False, - learn=True): - model = instance.to_model() - self.parent_solver.set_instance(model) - self.cplex = self.parent_solver._solver_model - x = self._get_features(instance) - - if self.ws_predictor is not None: - self.cplex.MIP_starts.delete() - ws = self.ws_predictor.predict(x) - if ws is not None: - _add_warm_start(self.cplex, ws) - - self.parent_solver.solve(tee=tee) + self.x_train = {} + self.y_train = {} + self.ws_predictors = {} - solution = np.array(self.cplex.solution.get_values()) - y = np.transpose(np.vstack((solution, 1 - solution))) - self._update_training_set(x, y) - return y - - def transform(self, instance): + def solve(self, instance, tee=False): + # Convert instance into concrete model model = instance.to_model() - self.parent_solver.set_instance(model) - self.cplex = self.parent_solver._solver_model - return self._get_features(instance) - - def predict(self, instance): - pass - def _update_training_set(self, x, y): - if self.train_x is None: - self.train_x = x - self.train_y = y - else: - self.train_x = np.vstack((self.train_x, x)) - self.train_y = np.vstack((self.train_y, y)) - - def fit(self): - if self.ws_predictor is not None: - self.ws_predictor.fit(self.train_x, self.train_y) + # Split decision variables according to their category + transformer = PerVariableTransformer() + var_split = transformer.split_variables(instance, model) + + # Build x_test and update x_train + x_test = {} + for category in var_split.keys(): + var_index_pairs = var_split[category] + x = transformer.transform_instance(instance, var_index_pairs) + x_test[category] = x + if category not in self.x_train.keys(): + self.x_train[category] = x + else: + self.x_train[category] = np.vstack([self.x_train[category], x]) + + # Predict warm start + for category in var_split.keys(): + if category in self.ws_predictors.keys(): + var_index_pairs = var_split[category] + ws = self.ws_predictors[category].predict(x_test[category]) + assert ws.shape == (len(var_index_pairs), 2) + for i in range(len(var_index_pairs)): + var, index = var_index_pairs[i] + if ws[i,0] == 1: + var[index].value = 1 + elif ws[i,1] == 1: + var[index].value = 0 + + # Solve MILP + self.parent_solver.solve(model, tee=tee, warmstart=True) + + # Update y_train + for category in var_split.keys(): + var_index_pairs = var_split[category] + y = transformer.transform_solution(var_index_pairs) + if category not in self.y_train.keys(): + self.y_train[category] = y + else: + self.y_train[category] = np.vstack([self.y_train[category], y]) -def _add_warm_start(cplex, ws): - assert isinstance(ws, np.ndarray) - assert ws.shape == (cplex.variables.get_num(),) - indices, values = [], [] - for k in range(len(ws)): - if isfinite(ws[k]): - indices += [k] - values += [ws[k]] - print("Adding warm start with %d values" % len(indices)) - cplex.MIP_starts.add([indices, values], cplex.MIP_starts.effort_level.solve_MIP) + def fit(self, x_train_dict=None, y_train_dict=None): + if x_train_dict is None: + x_train_dict = self.x_train + y_train_dict = self.y_train + for category in x_train_dict.keys(): + x_train = x_train_dict[category] + y_train = y_train_dict[category] + self.ws_predictors[category] = WarmStartPredictor() + self.ws_predictors[category].fit(x_train, y_train) + def _solve(self, tee): + self.parent_solver.solve(tee=tee) \ No newline at end of file diff --git a/miplearn/tests/test_solver.py b/miplearn/tests/test_solver.py new file mode 100644 index 0000000..deea525 --- /dev/null +++ b/miplearn/tests/test_solver.py @@ -0,0 +1,16 @@ +# MIPLearn: A Machine-Learning Framework for Mixed-Integer Optimization +# Copyright (C) 2019-2020 Argonne National Laboratory. All rights reserved. +# Written by Alinson S. Xavier + +from miplearn import LearningSolver +from miplearn.problems.knapsack import KnapsackInstance2 + + +def test_solver(): + instance = KnapsackInstance2(weights=[23., 26., 20., 18.], + prices=[505., 352., 458., 220.], + capacity=67.) + solver = LearningSolver() + solver.solve(instance) + solver.fit() + solver.solve(instance) diff --git a/miplearn/tests/test_transformer.py b/miplearn/tests/test_transformer.py index 83d9171..a027dbd 100644 --- a/miplearn/tests/test_transformer.py +++ b/miplearn/tests/test_transformer.py @@ -2,19 +2,22 @@ # Copyright (C) 2019-2020 Argonne National Laboratory. All rights reserved. # Written by Alinson S. Xavier -from miplearn import Instance, LearningSolver from miplearn.transformers import PerVariableTransformer from miplearn.problems.knapsack import KnapsackInstance, KnapsackInstance2 import numpy as np import pyomo.environ as pe + def test_transform(): transformer = PerVariableTransformer() instance = KnapsackInstance(weights=[23., 26., 20., 18.], prices=[505., 352., 458., 220.], capacity=67.) model = instance.to_model() - + solver = pe.SolverFactory('cbc') + solver.options["threads"] = 1 + solver.solve(model) + var_split = transformer.split_variables(instance, model) var_split_expected = { "default": [ @@ -26,7 +29,7 @@ def test_transform(): } assert var_split == var_split_expected var_index_pairs = [(model.x, i) for i in range(4)] - + x_actual = transformer.transform_instance(instance, var_index_pairs) x_expected = np.array([ [67., 21.75, 23., 505.], @@ -34,24 +37,29 @@ def test_transform(): [67., 21.75, 20., 458.], [67., 21.75, 18., 220.], ]) - assert x_expected.tolist() == x_actual.tolist() - - solver = pe.SolverFactory('cplex') - solver.options["threads"] = 1 + assert x_expected.tolist() == np.round(x_actual, decimals=2).tolist() + solver.solve(model) - y_actual = transformer.transform_solution(var_index_pairs) - y_expected = np.array([1., 0., 1., 1.]) + y_expected = np.array([ + [0., 1.], + [1., 0.], + [0., 1.], + [0., 1.], + ]) assert y_actual.tolist() == y_expected.tolist() - - + + def test_transform_with_categories(): transformer = PerVariableTransformer() instance = KnapsackInstance2(weights=[23., 26., 20., 18.], prices=[505., 352., 458., 220.], capacity=67.) model = instance.to_model() - + solver = pe.SolverFactory('cbc') + solver.options["threads"] = 1 + solver.solve(model) + var_split = transformer.split_variables(instance, model) var_split_expected = { 0: [(model.x, 0)], @@ -63,13 +71,13 @@ def test_transform_with_categories(): var_index_pairs = var_split[0] x_actual = transformer.transform_instance(instance, var_index_pairs) - x_expected = np.array([[23., 26., 20., 18., 505., 352., 458., 220.]]) - assert x_expected.tolist() == x_actual.tolist() + x_expected = np.array([ + [23., 26., 20., 18., 505., 352., 458., 220.] + ]) + assert x_expected.tolist() == np.round(x_actual, decimals=2).tolist() - solver = pe.SolverFactory('cplex') - solver.options["threads"] = 1 solver.solve(model) - + y_actual = transformer.transform_solution(var_index_pairs) - y_expected = np.array([1.]) - assert y_actual.tolist() == y_expected.tolist() \ No newline at end of file + y_expected = np.array([[0., 1.]]) + assert y_actual.tolist() == y_expected.tolist() diff --git a/miplearn/transformers.py b/miplearn/transformers.py index 733eef4..e75f9a1 100644 --- a/miplearn/transformers.py +++ b/miplearn/transformers.py @@ -5,14 +5,16 @@ import numpy as np from pyomo.core import Var + class PerVariableTransformer: """ Class that converts a miplearn.Instance into a matrix of features that is suitable for training machine learning models that make one decision per decision variable. """ + def __init__(self): pass - + def transform_instance(self, instance, var_index_pairs): instance_features = self._get_instance_features(instance) variable_features = self._get_variable_features(instance, var_index_pairs) @@ -20,13 +22,15 @@ class PerVariableTransformer: np.hstack([instance_features, vf]) for vf in variable_features ]) - - def _get_instance_features(self, instance): + + @staticmethod + def _get_instance_features(instance): features = instance.get_instance_features() assert isinstance(features, np.ndarray) return features - def _get_variable_features(self, instance, var_index_pairs): + @staticmethod + def _get_variable_features(instance, var_index_pairs): features = [] expected_shape = None for (var, index) in var_index_pairs: @@ -39,19 +43,21 @@ class PerVariableTransformer: assert vf.shape == expected_shape features += [vf] return np.array(features) - - def transform_solution(self, var_index_pairs): + + @staticmethod + def transform_solution(var_index_pairs): y = [] for (var, index) in var_index_pairs: - y += [var[index].value] + y += [[1 - var[index].value, var[index].value]] return np.array(y) - - def split_variables(self, instance, model): + + @staticmethod + def split_variables(instance, model): result = {} for var in model.component_objects(Var): for index in var: category = instance.get_variable_category(var, index) if category not in result.keys(): result[category] = [] - result[category] += [(var,index)] - return result \ No newline at end of file + result[category] += [(var, index)] + return result diff --git a/miplearn/warmstart.py b/miplearn/warmstart.py index 0ccb46d..7156919 100644 --- a/miplearn/warmstart.py +++ b/miplearn/warmstart.py @@ -2,28 +2,42 @@ # Copyright (C) 2019-2020 Argonne National Laboratory. All rights reserved. # Written by Alinson S. Xavier -import tensorflow as tf -import tensorflow.keras as keras -from tensorflow.keras.models import Sequential -from tensorflow.keras.layers import Dense, Dropout, Flatten, Activation import numpy as np +from sklearn.pipeline import make_pipeline +from sklearn.linear_model import LogisticRegression +from sklearn.preprocessing import StandardScaler + class WarmStartPredictor: - def __init__(self, model=None, threshold=0.80): - self.model = model - self.threshold = threshold - - def fit(self, train_x, train_y): - pass - - def predict(self, x): - if self.model is None: return None - assert isinstance(x, np.ndarray) - y = self.model.predict(x) - n_vars = y.shape[0] - ws = np.array([float("nan")] * n_vars) - ws[y[:,0] > self.threshold] = 1.0 - ws[y[:,1] > self.threshold] = 0.0 - return ws - - \ No newline at end of file + def __init__(self, + thr_fix_zero=0.05, + thr_fix_one=0.95, + thr_predict=0.95): + self.model = None + self.thr_predict = thr_predict + self.thr_fix_zero = thr_fix_zero + self.thr_fix_one = thr_fix_one + + def fit(self, x_train, y_train): + assert isinstance(x_train, np.ndarray) + assert isinstance(y_train, np.ndarray) + assert y_train.shape[1] == 2 + assert y_train.shape[0] == x_train.shape[0] + y_hat = np.average(y_train[:, 1]) + if y_hat < self.thr_fix_zero or y_hat > self.thr_fix_one: + self.model = int(y_hat) + else: + self.model = make_pipeline(StandardScaler(), LogisticRegression()) + self.model.fit(x_train, y_train[:, 1].astype(int)) + + def predict(self, x_test): + assert isinstance(x_test, np.ndarray) + if isinstance(self.model, int): + p_test = np.array([[1 - self.model, self.model] + for _ in range(x_test.shape[0])]) + else: + p_test = self.model.predict_proba(x_test) + p_test[p_test < self.thr_predict] = 0 + p_test[p_test > 0] = 1 + p_test = p_test.astype(int) + return p_test diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..048b09c --- /dev/null +++ b/requirements.txt @@ -0,0 +1,4 @@ +pyomo +numpy +pytest +sklearn \ No newline at end of file diff --git a/setup.py b/setup.py index 008182a..a64f10f 100644 --- a/setup.py +++ b/setup.py @@ -7,5 +7,5 @@ setup( author='Alinson S. Xavier', author_email='axavier@anl.gov', packages=['miplearn'], - install_requires=['pyomo'], + install_requires=['pyomo', 'numpy', 'sklearn'], ) \ No newline at end of file