diff --git a/docs/customization/index.html b/docs/customization/index.html
index 67b3a52..b509b90 100644
--- a/docs/customization/index.html
+++ b/docs/customization/index.html
@@ -142,6 +142,8 @@
                 
             <li class="second-level"><a href="#selecting-solver-components">Selecting solver components</a></li>
                 
+            <li class="second-level"><a href="#adjusting-component-aggresiveness">Adjusting component aggresiveness</a></li>
+                
     </ul>
 </div></div>
         <div class="col-md-9" role="main">
@@ -183,6 +185,19 @@ solver = LearningSolver()
 
 # Replace the default LazyConstraintComponent by one with custom parameters 
 solver.add(LazyConstraintComponent(...))
+</code></pre>
+
+<h2 id="adjusting-component-aggresiveness">Adjusting component aggresiveness</h2>
+<p>The aggressiveness of classification components (such as <code>PrimalSolutionComponent</code> and <code>LazyConstraintComponent</code>) can
+be adjusted through the <code>threshold</code> constructor argument. Internally, these components ask the ML models how confident
+they are on each prediction (through the <code>predict_proba</code> method in the sklearn API), and only take into account
+predictions which have probabilities above the threshold. Lowering a component's threshold increases its aggresiveness,
+while raising a component's threshold makes it more conservative. </p>
+<p>MIPLearn also includes <code>MinPrecisionThreshold</code>, a dynamic threshold which adjusts itself automatically during training
+to achieve a minimum desired true positive rate (also known as precision). The example below shows how to initialize
+a <code>PrimalSolutionComponent</code> which achieves 95% precision, possibly at the cost of a lower recall. To make the component
+more aggressive, this precision may be lowered.</p>
+<pre><code class="python">comp = PrimalSolutionComponent(threshold=MinPrecisionThreshold(0.98))
 </code></pre></div>
         
         
diff --git a/docs/index.html b/docs/index.html
index 01d8094..6cdab3c 100644
--- a/docs/index.html
+++ b/docs/index.html
@@ -273,5 +273,5 @@
 
 <!--
 MkDocs version : 1.1
-Build Date UTC : 2020-04-14 13:54:25
+Build Date UTC : 2020-05-05 18:08:32
 -->
diff --git a/docs/search/search_index.json b/docs/search/search_index.json
index 7e09197..dcb1217 100644
--- a/docs/search/search_index.json
+++ b/docs/search/search_index.json
@@ -1 +1 @@
-{"config":{"lang":["en"],"prebuild_index":false,"separator":"[\\s\\-]+"},"docs":[{"location":"","text":"MIPLearn MIPLearn is an extensible framework for Learning-Enhanced Mixed-Integer Optimization , an approach targeted at discrete optimization problems that need to be repeatedly solved with only minor changes to input data. The package uses Machine Learning (ML) to automatically identify patterns in previously solved instances of the problem, or in the solution process itself, and produces hints that can guide a conventional MIP solver towards the optimal solution faster. For particular classes of problems, this approach has been shown to provide significant performance benefits (see benchmark results and references for more details). Features MIPLearn proposes a flexible problem specification format, which allows users to describe their particular optimization problems to a Learning-Enhanced MIP solver, both from the MIP perspective and from the ML perspective, without making any assumptions on the problem being modeled, the mathematical formulation of the problem, or ML encoding. While the format is very flexible, some constraints are enforced to ensure that it is usable by an actual solver. MIPLearn provides a reference implementation of a Learning-Enhanced Solver , which can use the above problem specification format to automatically predict, based on previously solved instances, a number of hints to accelerate MIP performance. Currently, the reference solver is able to predict: (i) partial solutions which are likely to work well as MIP starts; (ii) an initial set of lazy constraints to enforce; (iii) variable branching priorities to accelerate the exploration of the branch-and-bound tree; (iv) the optimal objective value based on the solution to the LP relaxation. The usage of the solver is very straightforward. The most suitable ML models are automatically selected, trained, cross-validated and applied to the problem with no user intervention. MIPLearn provides a set of benchmark problems and random instance generators, covering applications from different domains, which can be used to quickly evaluate new learning-enhanced MIP techniques in a measurable and reproducible way. MIPLearn is customizable and extensible . For MIP and ML researchers exploring new techniques to accelerate MIP performance based on historical data, each component of the reference solver can be individually replaced, extended or customized. Documentation Installation and typical usage Benchmark utilities Benchmark problems, challenges and results Customizing the solver License, authors, references and acknowledgements Souce Code https://github.com/ANL-CEEESA/MIPLearn","title":"Home"},{"location":"#miplearn","text":"MIPLearn is an extensible framework for Learning-Enhanced Mixed-Integer Optimization , an approach targeted at discrete optimization problems that need to be repeatedly solved with only minor changes to input data. The package uses Machine Learning (ML) to automatically identify patterns in previously solved instances of the problem, or in the solution process itself, and produces hints that can guide a conventional MIP solver towards the optimal solution faster. For particular classes of problems, this approach has been shown to provide significant performance benefits (see benchmark results and references for more details).","title":"MIPLearn"},{"location":"#features","text":"MIPLearn proposes a flexible problem specification format, which allows users to describe their particular optimization problems to a Learning-Enhanced MIP solver, both from the MIP perspective and from the ML perspective, without making any assumptions on the problem being modeled, the mathematical formulation of the problem, or ML encoding. While the format is very flexible, some constraints are enforced to ensure that it is usable by an actual solver. MIPLearn provides a reference implementation of a Learning-Enhanced Solver , which can use the above problem specification format to automatically predict, based on previously solved instances, a number of hints to accelerate MIP performance. Currently, the reference solver is able to predict: (i) partial solutions which are likely to work well as MIP starts; (ii) an initial set of lazy constraints to enforce; (iii) variable branching priorities to accelerate the exploration of the branch-and-bound tree; (iv) the optimal objective value based on the solution to the LP relaxation. The usage of the solver is very straightforward. The most suitable ML models are automatically selected, trained, cross-validated and applied to the problem with no user intervention. MIPLearn provides a set of benchmark problems and random instance generators, covering applications from different domains, which can be used to quickly evaluate new learning-enhanced MIP techniques in a measurable and reproducible way. MIPLearn is customizable and extensible . For MIP and ML researchers exploring new techniques to accelerate MIP performance based on historical data, each component of the reference solver can be individually replaced, extended or customized.","title":"Features"},{"location":"#documentation","text":"Installation and typical usage Benchmark utilities Benchmark problems, challenges and results Customizing the solver License, authors, references and acknowledgements","title":"Documentation"},{"location":"#souce-code","text":"https://github.com/ANL-CEEESA/MIPLearn","title":"Souce Code"},{"location":"about/","text":"About Authors Alinson S. Xavier, Argonne National Laboratory < axavier@anl.gov > Feng Qiu, Argonne National Laboratory < fqiu@anl.gov > Acknowledgements Based upon work supported by Laboratory Directed Research and Development (LDRD) funding from Argonne National Laboratory, provided by the Director, Office of Science, of the U.S. Department of Energy under Contract No. DE-AC02-06CH11357. References Learning to Solve Large-Scale Security-Constrained Unit Commitment Problems. Alinson S. Xavier, Feng Qiu, Shabbir Ahmed . INFORMS Journal on Computing (to appear). ArXiv:1902:01696 License MIPLearn, an extensible framework for Learning-Enhanced Mixed-Integer Optimization Copyright \u00a9 2020, UChicago Argonne, LLC. All Rights Reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \"AS IS\" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.","title":"About"},{"location":"about/#about","text":"","title":"About"},{"location":"about/#authors","text":"Alinson S. Xavier, Argonne National Laboratory < axavier@anl.gov > Feng Qiu, Argonne National Laboratory < fqiu@anl.gov >","title":"Authors"},{"location":"about/#acknowledgements","text":"Based upon work supported by Laboratory Directed Research and Development (LDRD) funding from Argonne National Laboratory, provided by the Director, Office of Science, of the U.S. Department of Energy under Contract No. DE-AC02-06CH11357.","title":"Acknowledgements"},{"location":"about/#references","text":"Learning to Solve Large-Scale Security-Constrained Unit Commitment Problems. Alinson S. Xavier, Feng Qiu, Shabbir Ahmed . INFORMS Journal on Computing (to appear). ArXiv:1902:01696","title":"References"},{"location":"about/#license","text":"MIPLearn, an extensible framework for Learning-Enhanced Mixed-Integer Optimization Copyright \u00a9 2020, UChicago Argonne, LLC. All Rights Reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \"AS IS\" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.","title":"License"},{"location":"benchmark/","text":"Benchmarks Utilities Using BenchmarkRunner MIPLearn provides the utility class BenchmarkRunner , which simplifies the task of comparing the performance of different solvers. The snippet below shows its basic usage: from miplearn import BenchmarkRunner, LearningSolver # Create train and test instances train_instances = [...] test_instances = [...] # Training phase... training_solver = LearningSolver(...) training_solver.parallel_solve(train_instances, n_jobs=10) # Test phase... test_solvers = { \"Baseline\": LearningSolver(...), # each solver may have different parameters \"Strategy A\": LearningSolver(...), \"Strategy B\": LearningSolver(...), \"Strategy C\": LearningSolver(...), } benchmark = BenchmarkRunner(test_solvers) benchmark.fit(train_instances) benchmark.parallel_solve(test_instances, n_jobs=2) print(benchmark.raw_results()) The method fit trains the ML models for each individual solver. The method parallel_solve solves the test instances in parallel, and collects solver statistics such as running time and optimal value. Finally, raw_results produces a table of results (Pandas DataFrame) with the following columns: Solver, the name of the solver. Instance, the sequence number identifying the instance. Wallclock Time, the wallclock running time (in seconds) spent by the solver; Lower Bound, the best lower bound obtained by the solver; Upper Bound, the best upper bound obtained by the solver; Gap, the relative MIP integrality gap at the end of the optimization; Nodes, the number of explored branch-and-bound nodes. In addition to the above, there is also a \"Relative\" version of most columns, where the raw number is compared to the solver which provided the best performance. The Relative Wallclock Time for example, indicates how many times slower this run was when compared to the best time achieved by any solver when processing this instance. For example, if this run took 10 seconds, but the fastest solver took only 5 seconds to solve the same instance, the relative wallclock time would be 2. Saving and loading benchmark results When iteratively exploring new formulations, encoding and solver parameters, it is often desirable to avoid repeating parts of the benchmark suite. For example, if the baseline solver has not been changed, there is no need to evaluate its performance again and again when making small changes to the remaining solvers. BenchmarkRunner provides the methods save_results and load_results , which can be used to avoid this repetition, as the next example shows: # Benchmark baseline solvers and save results to a file. benchmark = BenchmarkRunner(baseline_solvers) benchmark.parallel_solve(test_instances) benchmark.save_results(\"baseline_results.csv\") # Benchmark remaining solvers, loading baseline results from file. benchmark = BenchmarkRunner(alternative_solvers) benchmark.load_results(\"baseline_results.csv\") benchmark.fit(training_instances) benchmark.parallel_solve(test_instances)","title":"Benchmark"},{"location":"benchmark/#benchmarks-utilities","text":"","title":"Benchmarks Utilities"},{"location":"benchmark/#using-benchmarkrunner","text":"MIPLearn provides the utility class BenchmarkRunner , which simplifies the task of comparing the performance of different solvers. The snippet below shows its basic usage: from miplearn import BenchmarkRunner, LearningSolver # Create train and test instances train_instances = [...] test_instances = [...] # Training phase... training_solver = LearningSolver(...) training_solver.parallel_solve(train_instances, n_jobs=10) # Test phase... test_solvers = { \"Baseline\": LearningSolver(...), # each solver may have different parameters \"Strategy A\": LearningSolver(...), \"Strategy B\": LearningSolver(...), \"Strategy C\": LearningSolver(...), } benchmark = BenchmarkRunner(test_solvers) benchmark.fit(train_instances) benchmark.parallel_solve(test_instances, n_jobs=2) print(benchmark.raw_results()) The method fit trains the ML models for each individual solver. The method parallel_solve solves the test instances in parallel, and collects solver statistics such as running time and optimal value. Finally, raw_results produces a table of results (Pandas DataFrame) with the following columns: Solver, the name of the solver. Instance, the sequence number identifying the instance. Wallclock Time, the wallclock running time (in seconds) spent by the solver; Lower Bound, the best lower bound obtained by the solver; Upper Bound, the best upper bound obtained by the solver; Gap, the relative MIP integrality gap at the end of the optimization; Nodes, the number of explored branch-and-bound nodes. In addition to the above, there is also a \"Relative\" version of most columns, where the raw number is compared to the solver which provided the best performance. The Relative Wallclock Time for example, indicates how many times slower this run was when compared to the best time achieved by any solver when processing this instance. For example, if this run took 10 seconds, but the fastest solver took only 5 seconds to solve the same instance, the relative wallclock time would be 2.","title":"Using BenchmarkRunner"},{"location":"benchmark/#saving-and-loading-benchmark-results","text":"When iteratively exploring new formulations, encoding and solver parameters, it is often desirable to avoid repeating parts of the benchmark suite. For example, if the baseline solver has not been changed, there is no need to evaluate its performance again and again when making small changes to the remaining solvers. BenchmarkRunner provides the methods save_results and load_results , which can be used to avoid this repetition, as the next example shows: # Benchmark baseline solvers and save results to a file. benchmark = BenchmarkRunner(baseline_solvers) benchmark.parallel_solve(test_instances) benchmark.save_results(\"baseline_results.csv\") # Benchmark remaining solvers, loading baseline results from file. benchmark = BenchmarkRunner(alternative_solvers) benchmark.load_results(\"baseline_results.csv\") benchmark.fit(training_instances) benchmark.parallel_solve(test_instances)","title":"Saving and loading benchmark results"},{"location":"customization/","text":"Customization Selecting the internal MIP solver By default, LearningSolver uses Gurobi as its internal MIP solver. Another supported solver is IBM ILOG CPLEX . To switch between solvers, use the solver constructor argument, as shown below. It is also possible to specify a time limit (in seconds) and a relative MIP gap tolerance. from miplearn import LearningSolver solver = LearningSolver(solver=\"cplex\", time_limit=300, gap_tolerance=1e-3) Selecting solver components LearningSolver is composed by a number of individual machine-learning components, each targeting a different part of the solution process. Each component can be individually enabled, disabled or customized. The following components are enabled by default: LazyConstraintComponent : Predicts which lazy constraint to initially enforce. ObjectiveValueComponent : Predicts the optimal value of the optimization problem, given the optimal solution to the LP relaxation. PrimalSolutionComponent : Predicts optimal values for binary decision variables. In heuristic mode, this component fixes the variables to their predicted values. In exact mode, the predicted values are provided to the solver as a (partial) MIP start. The following components are also available, but not enabled by default: BranchPriorityComponent : Predicts good branch priorities for decision variables. To create a LearningSolver with a specific set of components, the components constructor argument may be used, as the next example shows: # Create a solver without any components solver1 = LearningSolver(components=[]) # Create a solver with only two components solver2 = LearningSolver(components=[ LazyConstraintComponent(...), PrimalSolutionComponent(...), ]) It is also possible to add components to an existing solver using the solver.add method, as shown below. If the solver already holds another component of that type, the new component will replace the previous one. # Create solver with default components solver = LearningSolver() # Replace the default LazyConstraintComponent by one with custom parameters solver.add(LazyConstraintComponent(...))","title":"Customization"},{"location":"customization/#customization","text":"","title":"Customization"},{"location":"customization/#selecting-the-internal-mip-solver","text":"By default, LearningSolver uses Gurobi as its internal MIP solver. Another supported solver is IBM ILOG CPLEX . To switch between solvers, use the solver constructor argument, as shown below. It is also possible to specify a time limit (in seconds) and a relative MIP gap tolerance. from miplearn import LearningSolver solver = LearningSolver(solver=\"cplex\", time_limit=300, gap_tolerance=1e-3)","title":"Selecting the internal MIP solver"},{"location":"customization/#selecting-solver-components","text":"LearningSolver is composed by a number of individual machine-learning components, each targeting a different part of the solution process. Each component can be individually enabled, disabled or customized. The following components are enabled by default: LazyConstraintComponent : Predicts which lazy constraint to initially enforce. ObjectiveValueComponent : Predicts the optimal value of the optimization problem, given the optimal solution to the LP relaxation. PrimalSolutionComponent : Predicts optimal values for binary decision variables. In heuristic mode, this component fixes the variables to their predicted values. In exact mode, the predicted values are provided to the solver as a (partial) MIP start. The following components are also available, but not enabled by default: BranchPriorityComponent : Predicts good branch priorities for decision variables. To create a LearningSolver with a specific set of components, the components constructor argument may be used, as the next example shows: # Create a solver without any components solver1 = LearningSolver(components=[]) # Create a solver with only two components solver2 = LearningSolver(components=[ LazyConstraintComponent(...), PrimalSolutionComponent(...), ]) It is also possible to add components to an existing solver using the solver.add method, as shown below. If the solver already holds another component of that type, the new component will replace the previous one. # Create solver with default components solver = LearningSolver() # Replace the default LazyConstraintComponent by one with custom parameters solver.add(LazyConstraintComponent(...))","title":"Selecting solver components"},{"location":"problems/","text":"Benchmark Problems, Challenges and Results MIPLearn provides a selection of benchmark problems and random instance generators, covering applications from different fields, that can be used to evaluate new learning-enhanced MIP techniques in a measurable and reproducible way. In this page, we describe these problems, the included instance generators, and we present some benchmark results for LearningSolver with default parameters. Preliminaries Benchmark challenges When evaluating the performance of a conventional MIP solver, benchmark sets , such as MIPLIB and TSPLIB, are typically used. The performance of newly proposed solvers or solution techniques are typically measured as the average (or total) running time the solver takes to solve the entire benchmark set. For Learning-Enhanced MIP solvers, it is also necessary to specify what instances should the solver be trained on (the training instances ) before solving the actual set of instances we are interested in (the test instances ). If the training instances are very similar to the test instances, we would expect a Learning-Enhanced Solver to present stronger perfomance benefits. In MIPLearn, each optimization problem comes with a set of benchmark challenges , which specify how should the training and test instances be generated. The first challenges are typically easier, in the sense that training and test instances are very similar. Later challenges gradually make the sets more distinct, and therefore harder to learn from. Baseline results To illustrate the performance of LearningSolver , and to set a baseline for newly proposed techniques, we present in this page, for each benchmark challenge, a small set of computational results measuring the solution speed of the solver and the solution quality with default parameters. For more detailed computational studies, see references . We compare three solvers: baseline: Gurobi 9.0 with default settings (a conventional state-of-the-art MIP solver) ml-exact: LearningSolver with default settings, using Gurobi 9.0 as internal MIP solver ml-heuristic: Same as above, but with mode=\"heuristic\" All experiments presented here were performed on a Linux server (Ubuntu Linux 18.04 LTS) with Intel Xeon Gold 6230s (2 processors, 40 cores, 80 threads) and 256 GB RAM (DDR4, 2933 MHz). All solvers were restricted to use 4 threads, with no time limits, and 10 instances were solved simultaneously at a time. Maximum Weight Stable Set Problem Problem definition Given a simple undirected graph $G=(V,E)$ and weights $w \\in \\mathbb{R}^V$, the problem is to find a stable set $S \\subseteq V$ that maximizes $ \\sum_{v \\in V} w_v$. We recall that a subset $S \\subseteq V$ is a stable set if no two vertices of $S$ are adjacent. This is one of Karp's 21 NP-complete problems. Random instance generator The class MaxWeightStableSetGenerator can be used to generate random instances of this problem, with user-specified probability distributions. When the constructor parameter fix_graph=True is provided, one random Erd\u0151s-R\u00e9nyi graph $G_{n,p}$ is generated during the constructor, where $n$ and $p$ are sampled from user-provided probability distributions n and p . To generate each instance, the generator independently samples each $w_v$ from the user-provided probability distribution w . When fix_graph=False , a new random graph is generated for each instance, while the remaining parameters are sampled in the same way. Challenge A Fixed random Erd\u0151s-R\u00e9nyi graph $G_{n,p}$ with $n=200$ and $p=5\\%$ Random vertex weights $w_v \\sim U(100, 150)$ 500 training instances, 50 test instances MaxWeightStableSetGenerator(w=uniform(loc=100., scale=50.), n=randint(low=200, high=201), p=uniform(loc=0.05, scale=0.0), fix_graph=True) Traveling Salesman Problem Problem definition Given a list of cities and the distance between each pair of cities, the problem asks for the shortest route starting at the first city, visiting each other city exactly once, then returning to the first city. This problem is a generalization of the Hamiltonian path problem, one of Karp's 21 NP-complete problems. Random problem generator The class TravelingSalesmanGenerator can be used to generate random instances of this problem. Initially, the generator creates $n$ cities $(x_1,y_1),\\ldots,(x_n,y_n) \\in \\mathbb{R}^2$, where $n, x_i$ and $y_i$ are sampled independently from the provided probability distributions n , x and y . For each pair of cities $(i,j)$, the distance $d_{i,j}$ between them is set to: d_{i,j} = \\gamma_{i,j} \\sqrt{(x_i-x_j)^2 + (y_i - y_j)^2} where $\\gamma_{i,j}$ is sampled from the distribution gamma . If fix_cities=True is provided, 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. Challenge A Fixed list of 350 cities in the $[0, 1000]^2$ square $\\gamma_{i,j} \\sim U(0.95, 1.05)$ 500 training instances, 50 test instances 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, round=True, ) Multidimensional 0-1 Knapsack Problem Problem definition Given a set of $n$ items and $m$ types of resources (also called knapsacks ), the problem is to find a subset of items that maximizes profit without consuming more resources than it is available. More precisely, the problem is: \\begin{align*} \\text{maximize} & \\sum_{j=1}^n p_j x_j \\\\ \\text{subject to} & \\sum_{j=1}^n w_{ij} x_j \\leq b_i & \\forall i=1,\\ldots,m \\\\ & x_j \\in \\{0,1\\} & \\forall j=1,\\ldots,n \\end{align*} Random instance generator The class MultiKnapsackGenerator can be used to generate random instances of this problem. The number of items $n$ and knapsacks $m$ are sampled from the user-provided probability distributions n and m . The weights $w_{ij}$ are sampled independently from the provided distribution w . The capacity of knapsack $i$ is set to b_i = \\alpha_i \\sum_{j=1}^n w_{ij} where $\\alpha_i$, the tightness ratio, is sampled from the provided probability distribution alpha . To make the instances more challenging, the costs of the items are linearly correlated to their average weights. More specifically, the price of each item $j$ is set to: p_j = \\sum_{i=1}^m \\frac{w_{ij}}{m} + K u_j, where $K$, the correlation coefficient, and $u_j$, the correlation multiplier, are sampled from the provided probability distributions K and u . If fix_w=True is provided, then $w_{ij}$ are kept the same in all generated instances. This also implies that $n$ and $m$ are kept fixed. Although the prices and capacities are derived from $w_{ij}$, as long as u and K are not constants, the generated instances will still not be completely identical. If a probability distribution w_jitter is provided, then item weights will be set to $w_{ij} \\gamma_{ij}$ where $\\gamma_{ij}$ is sampled from w_jitter . When combined with fix_w=True , this argument may be used to generate instances where the weight of each item is roughly the same, but not exactly identical, across all instances. The prices of the items and the capacities of the knapsacks will be calculated as above, but using these perturbed weights instead. By default, all generated prices, weights and capacities are rounded to the nearest integer number. If round=False is provided, this rounding will be disabled. References Freville, Arnaud, and G\u00e9rard Plateau. An efficient preprocessing procedure for the multidimensional 0\u20131 knapsack problem. Discrete applied mathematics 49.1-3 (1994): 189-212. Fr\u00e9ville, Arnaud. The multidimensional 0\u20131 knapsack problem: An overview. European Journal of Operational Research 155.1 (2004): 1-21. Challenge A 250 variables, 10 constraints, fixed weights $w \\sim U(0, 1000), \\gamma \\sim U(0.95, 1.05)$ $K = 500, u \\sim U(0, 1), \\alpha = 0.25$ 500 training instances, 50 test instances MultiKnapsackGenerator(n=randint(low=250, high=251), m=randint(low=10, high=11), w=uniform(loc=0.0, scale=1000.0), K=uniform(loc=500.0, scale=0.0), u=uniform(loc=0.0, scale=1.0), alpha=uniform(loc=0.25, scale=0.0), fix_w=True, w_jitter=uniform(loc=0.95, scale=0.1), )","title":"Problems"},{"location":"problems/#benchmark-problems-challenges-and-results","text":"MIPLearn provides a selection of benchmark problems and random instance generators, covering applications from different fields, that can be used to evaluate new learning-enhanced MIP techniques in a measurable and reproducible way. In this page, we describe these problems, the included instance generators, and we present some benchmark results for LearningSolver with default parameters.","title":"Benchmark Problems, Challenges and Results"},{"location":"problems/#preliminaries","text":"","title":"Preliminaries"},{"location":"problems/#benchmark-challenges","text":"When evaluating the performance of a conventional MIP solver, benchmark sets , such as MIPLIB and TSPLIB, are typically used. The performance of newly proposed solvers or solution techniques are typically measured as the average (or total) running time the solver takes to solve the entire benchmark set. For Learning-Enhanced MIP solvers, it is also necessary to specify what instances should the solver be trained on (the training instances ) before solving the actual set of instances we are interested in (the test instances ). If the training instances are very similar to the test instances, we would expect a Learning-Enhanced Solver to present stronger perfomance benefits. In MIPLearn, each optimization problem comes with a set of benchmark challenges , which specify how should the training and test instances be generated. The first challenges are typically easier, in the sense that training and test instances are very similar. Later challenges gradually make the sets more distinct, and therefore harder to learn from.","title":"Benchmark challenges"},{"location":"problems/#baseline-results","text":"To illustrate the performance of LearningSolver , and to set a baseline for newly proposed techniques, we present in this page, for each benchmark challenge, a small set of computational results measuring the solution speed of the solver and the solution quality with default parameters. For more detailed computational studies, see references . We compare three solvers: baseline: Gurobi 9.0 with default settings (a conventional state-of-the-art MIP solver) ml-exact: LearningSolver with default settings, using Gurobi 9.0 as internal MIP solver ml-heuristic: Same as above, but with mode=\"heuristic\" All experiments presented here were performed on a Linux server (Ubuntu Linux 18.04 LTS) with Intel Xeon Gold 6230s (2 processors, 40 cores, 80 threads) and 256 GB RAM (DDR4, 2933 MHz). All solvers were restricted to use 4 threads, with no time limits, and 10 instances were solved simultaneously at a time.","title":"Baseline results"},{"location":"problems/#maximum-weight-stable-set-problem","text":"","title":"Maximum Weight Stable Set Problem"},{"location":"problems/#problem-definition","text":"Given a simple undirected graph $G=(V,E)$ and weights $w \\in \\mathbb{R}^V$, the problem is to find a stable set $S \\subseteq V$ that maximizes $ \\sum_{v \\in V} w_v$. We recall that a subset $S \\subseteq V$ is a stable set if no two vertices of $S$ are adjacent. This is one of Karp's 21 NP-complete problems.","title":"Problem definition"},{"location":"problems/#random-instance-generator","text":"The class MaxWeightStableSetGenerator can be used to generate random instances of this problem, with user-specified probability distributions. When the constructor parameter fix_graph=True is provided, one random Erd\u0151s-R\u00e9nyi graph $G_{n,p}$ is generated during the constructor, where $n$ and $p$ are sampled from user-provided probability distributions n and p . To generate each instance, the generator independently samples each $w_v$ from the user-provided probability distribution w . When fix_graph=False , a new random graph is generated for each instance, while the remaining parameters are sampled in the same way.","title":"Random instance generator"},{"location":"problems/#challenge-a","text":"Fixed random Erd\u0151s-R\u00e9nyi graph $G_{n,p}$ with $n=200$ and $p=5\\%$ Random vertex weights $w_v \\sim U(100, 150)$ 500 training instances, 50 test instances MaxWeightStableSetGenerator(w=uniform(loc=100., scale=50.), n=randint(low=200, high=201), p=uniform(loc=0.05, scale=0.0), fix_graph=True)","title":"Challenge A"},{"location":"problems/#traveling-salesman-problem","text":"","title":"Traveling Salesman Problem"},{"location":"problems/#problem-definition_1","text":"Given a list of cities and the distance between each pair of cities, the problem asks for the shortest route starting at the first city, visiting each other city exactly once, then returning to the first city. This problem is a generalization of the Hamiltonian path problem, one of Karp's 21 NP-complete problems.","title":"Problem definition"},{"location":"problems/#random-problem-generator","text":"The class TravelingSalesmanGenerator can be used to generate random instances of this problem. Initially, the generator creates $n$ cities $(x_1,y_1),\\ldots,(x_n,y_n) \\in \\mathbb{R}^2$, where $n, x_i$ and $y_i$ are sampled independently from the provided probability distributions n , x and y . For each pair of cities $(i,j)$, the distance $d_{i,j}$ between them is set to: d_{i,j} = \\gamma_{i,j} \\sqrt{(x_i-x_j)^2 + (y_i - y_j)^2} where $\\gamma_{i,j}$ is sampled from the distribution gamma . If fix_cities=True is provided, 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.","title":"Random problem generator"},{"location":"problems/#challenge-a_1","text":"Fixed list of 350 cities in the $[0, 1000]^2$ square $\\gamma_{i,j} \\sim U(0.95, 1.05)$ 500 training instances, 50 test instances 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, round=True, )","title":"Challenge A"},{"location":"problems/#multidimensional-0-1-knapsack-problem","text":"","title":"Multidimensional 0-1 Knapsack Problem"},{"location":"problems/#problem-definition_2","text":"Given a set of $n$ items and $m$ types of resources (also called knapsacks ), the problem is to find a subset of items that maximizes profit without consuming more resources than it is available. More precisely, the problem is: \\begin{align*} \\text{maximize} & \\sum_{j=1}^n p_j x_j \\\\ \\text{subject to} & \\sum_{j=1}^n w_{ij} x_j \\leq b_i & \\forall i=1,\\ldots,m \\\\ & x_j \\in \\{0,1\\} & \\forall j=1,\\ldots,n \\end{align*}","title":"Problem definition"},{"location":"problems/#random-instance-generator_1","text":"The class MultiKnapsackGenerator can be used to generate random instances of this problem. The number of items $n$ and knapsacks $m$ are sampled from the user-provided probability distributions n and m . The weights $w_{ij}$ are sampled independently from the provided distribution w . The capacity of knapsack $i$ is set to b_i = \\alpha_i \\sum_{j=1}^n w_{ij} where $\\alpha_i$, the tightness ratio, is sampled from the provided probability distribution alpha . To make the instances more challenging, the costs of the items are linearly correlated to their average weights. More specifically, the price of each item $j$ is set to: p_j = \\sum_{i=1}^m \\frac{w_{ij}}{m} + K u_j, where $K$, the correlation coefficient, and $u_j$, the correlation multiplier, are sampled from the provided probability distributions K and u . If fix_w=True is provided, then $w_{ij}$ are kept the same in all generated instances. This also implies that $n$ and $m$ are kept fixed. Although the prices and capacities are derived from $w_{ij}$, as long as u and K are not constants, the generated instances will still not be completely identical. If a probability distribution w_jitter is provided, then item weights will be set to $w_{ij} \\gamma_{ij}$ where $\\gamma_{ij}$ is sampled from w_jitter . When combined with fix_w=True , this argument may be used to generate instances where the weight of each item is roughly the same, but not exactly identical, across all instances. The prices of the items and the capacities of the knapsacks will be calculated as above, but using these perturbed weights instead. By default, all generated prices, weights and capacities are rounded to the nearest integer number. If round=False is provided, this rounding will be disabled. References Freville, Arnaud, and G\u00e9rard Plateau. An efficient preprocessing procedure for the multidimensional 0\u20131 knapsack problem. Discrete applied mathematics 49.1-3 (1994): 189-212. Fr\u00e9ville, Arnaud. The multidimensional 0\u20131 knapsack problem: An overview. European Journal of Operational Research 155.1 (2004): 1-21.","title":"Random instance generator"},{"location":"problems/#challenge-a_2","text":"250 variables, 10 constraints, fixed weights $w \\sim U(0, 1000), \\gamma \\sim U(0.95, 1.05)$ $K = 500, u \\sim U(0, 1), \\alpha = 0.25$ 500 training instances, 50 test instances MultiKnapsackGenerator(n=randint(low=250, high=251), m=randint(low=10, high=11), w=uniform(loc=0.0, scale=1000.0), K=uniform(loc=500.0, scale=0.0), u=uniform(loc=0.0, scale=1.0), alpha=uniform(loc=0.25, scale=0.0), fix_w=True, w_jitter=uniform(loc=0.95, scale=0.1), )","title":"Challenge A"},{"location":"usage/","text":"Usage Installation MIPLearn is mainly written in Python, with some components written in Julia. For this reason, both Python 3.6+ and Julia 1.3+ are required. A mixed-integer solver is also required, and its Python bindings must be properly installed. Supported solvers are CPLEX and Gurobi. Optimization problems currently need to be specified in the Pyomo modeling language. A JuMP interface to the package is currently under development. To install MIPLearn, run the following commands: git clone https://github.com/ANL-CEEESA/MIPLearn.git cd MIPLearn make install After installation, the package miplearn should become available to Python. It can be imported as follows: import miplearn Note To install MIPLearn in another Python environment, switch to that enviroment before running make install . To install the package in development mode, run make develop instead. Using LearningSolver The main class provided by this package is LearningSolver , a learning-enhanced MIP solver which uses information from previously solved instances to accelerate the solution of new instances. The following example shows its basic usage: from miplearn import LearningSolver # List of user-provided instances training_instances = [...] test_instances = [...] # Create solver solver = LearningSolver() # Solve all training instances for instance in training_instances: solver.solve(instance) # Learn from training instances solver.fit(training_instances) # Solve all test instances for instance in test_instances: solver.solve(instance) In this example, we have two lists of user-provided instances: training_instances and test_instances . We start by solving all training instances. Since there is no historical information available at this point, the instances will be processed from scratch, with no ML acceleration. After solving each instance, the solver stores within each instance object the optimal solution, the optimal objective value, and other information that can be used to accelerate future solves. After all training instances are solved, we call solver.fit(training_instances) . This instructs the solver to train all its internal machine-learning models based on the solutions of the (solved) trained instances. Subsequent calls to solver.solve(instance) will automatically use the trained Machine Learning models to accelerate the solution process. Describing problem instances Instances to be solved by LearningSolver must derive from the abstract class miplearn.Instance . The following three abstract methods must be implemented: instance.to_model() , which returns a concrete Pyomo model corresponding to the instance; instance.get_instance_features() , which returns a 1-dimensional Numpy array of (numerical) features describing the entire instance; instance.get_variable_features(var_name, index) , which returns a 1-dimensional array of (numerical) features describing a particular decision variable. The first method is used by LearningSolver to construct a concrete Pyomo model, which will be provided to the internal MIP solver. The second and third methods provide an encoding of the instance, which can be used by the ML models to make predictions. In the knapsack problem, for example, an implementation may decide to provide as instance features the average weights, average prices, number of items and the size of the knapsack. The weight and the price of each individual item could be provided as variable features. See src/python/miplearn/problems/knapsack.py for a concrete example. An optional method which can be implemented is instance.get_variable_category(var_name, index) , which returns a category (a string, an integer or any hashable type) for each decision variable. If two variables have the same category, LearningSolver will use the same internal ML model to predict the values of both variables. By default, all variables belong to the \"default\" category, and therefore only one ML model is used for all variables. If the returned category is None , ML predictors will ignore the variable. It is not necessary to have a one-to-one correspondence between features and problem instances. One important (and deliberate) limitation of MIPLearn, however, is that get_instance_features() must always return arrays of same length for all relevant instances of the problem. Similarly, get_variable_features(var_name, index) must also always return arrays of same length for all variables in each category. It is up to the user to decide how to encode variable-length characteristics of the problem into fixed-length vectors. In graph problems, for example, graph embeddings can be used to reduce the (variable-length) lists of nodes and edges into a fixed-length structure that still preserves some properties of the graph. Different instance encodings may have significant impact on performance. Obtaining heuristic solutions By default, LearningSolver uses Machine Learning to accelerate the MIP solution process, while maintaining all optimality guarantees provided by the MIP solver. In the default mode of operation, for example, predicted optimal solutions are used only as MIP starts. For more significant performance benefits, LearningSolver can also be configured to place additional trust in the Machine Learning predictors, by using the mode=\"heuristic\" constructor argument. When operating in this mode, if a ML model is statistically shown (through stratified k-fold cross validation ) to have exceptionally high accuracy, the solver may decide to restrict the search space based on its predictions. The parts of the solution which the ML models cannot predict accurately will still be explored using traditional (branch-and-bound) methods. For particular applications, this mode has been shown to quickly produce optimal or near-optimal solutions (see references and benchmark results ). Danger The heuristic mode provides no optimality guarantees, and therefore should only be used if the solver is first trained on a large and representative set of training instances. Training on a small or non-representative set of instances may produce low-quality solutions, or make the solver incorrectly classify new instances as infeasible. Saving and loading solver state After solving a large number of training instances, it may be desirable to save the current state of LearningSolver to disk, so that the solver can still use the acquired knowledge after the application restarts. This can be accomplished by using the standard pickle module, as the following example illustrates: from miplearn import LearningSolver import pickle # Solve training instances training_instances = [...] solver = LearningSolver() for instance in training_instances: solver.solve(instance) # Train machine-learning models solver.fit(training_instances) # Save trained solver to disk pickle.dump(solver, open(\"solver.pickle\", \"wb\")) # Application restarts... # Load trained solver from disk solver = pickle.load(open(\"solver.pickle\", \"rb\")) # Solve additional instances test_instances = [...] for instance in test_instances: solver.solve(instance) Solving training instances in parallel In many situations, training and test instances can be solved in parallel to accelerate the training process. LearningSolver provides the method parallel_solve(instances) to easily achieve this: from miplearn import LearningSolver training_instances = [...] solver = LearningSolver() solver.parallel_solve(training_instances, n_jobs=4) solver.fit(training_instances) # Test phase... test_instances = [...] solver.parallel_solve(test_instances) Current Limitations Only binary and continuous decision variables are currently supported.","title":"Usage"},{"location":"usage/#usage","text":"","title":"Usage"},{"location":"usage/#installation","text":"MIPLearn is mainly written in Python, with some components written in Julia. For this reason, both Python 3.6+ and Julia 1.3+ are required. A mixed-integer solver is also required, and its Python bindings must be properly installed. Supported solvers are CPLEX and Gurobi. Optimization problems currently need to be specified in the Pyomo modeling language. A JuMP interface to the package is currently under development. To install MIPLearn, run the following commands: git clone https://github.com/ANL-CEEESA/MIPLearn.git cd MIPLearn make install After installation, the package miplearn should become available to Python. It can be imported as follows: import miplearn Note To install MIPLearn in another Python environment, switch to that enviroment before running make install . To install the package in development mode, run make develop instead.","title":"Installation"},{"location":"usage/#using-learningsolver","text":"The main class provided by this package is LearningSolver , a learning-enhanced MIP solver which uses information from previously solved instances to accelerate the solution of new instances. The following example shows its basic usage: from miplearn import LearningSolver # List of user-provided instances training_instances = [...] test_instances = [...] # Create solver solver = LearningSolver() # Solve all training instances for instance in training_instances: solver.solve(instance) # Learn from training instances solver.fit(training_instances) # Solve all test instances for instance in test_instances: solver.solve(instance) In this example, we have two lists of user-provided instances: training_instances and test_instances . We start by solving all training instances. Since there is no historical information available at this point, the instances will be processed from scratch, with no ML acceleration. After solving each instance, the solver stores within each instance object the optimal solution, the optimal objective value, and other information that can be used to accelerate future solves. After all training instances are solved, we call solver.fit(training_instances) . This instructs the solver to train all its internal machine-learning models based on the solutions of the (solved) trained instances. Subsequent calls to solver.solve(instance) will automatically use the trained Machine Learning models to accelerate the solution process.","title":"Using LearningSolver"},{"location":"usage/#describing-problem-instances","text":"Instances to be solved by LearningSolver must derive from the abstract class miplearn.Instance . The following three abstract methods must be implemented: instance.to_model() , which returns a concrete Pyomo model corresponding to the instance; instance.get_instance_features() , which returns a 1-dimensional Numpy array of (numerical) features describing the entire instance; instance.get_variable_features(var_name, index) , which returns a 1-dimensional array of (numerical) features describing a particular decision variable. The first method is used by LearningSolver to construct a concrete Pyomo model, which will be provided to the internal MIP solver. The second and third methods provide an encoding of the instance, which can be used by the ML models to make predictions. In the knapsack problem, for example, an implementation may decide to provide as instance features the average weights, average prices, number of items and the size of the knapsack. The weight and the price of each individual item could be provided as variable features. See src/python/miplearn/problems/knapsack.py for a concrete example. An optional method which can be implemented is instance.get_variable_category(var_name, index) , which returns a category (a string, an integer or any hashable type) for each decision variable. If two variables have the same category, LearningSolver will use the same internal ML model to predict the values of both variables. By default, all variables belong to the \"default\" category, and therefore only one ML model is used for all variables. If the returned category is None , ML predictors will ignore the variable. It is not necessary to have a one-to-one correspondence between features and problem instances. One important (and deliberate) limitation of MIPLearn, however, is that get_instance_features() must always return arrays of same length for all relevant instances of the problem. Similarly, get_variable_features(var_name, index) must also always return arrays of same length for all variables in each category. It is up to the user to decide how to encode variable-length characteristics of the problem into fixed-length vectors. In graph problems, for example, graph embeddings can be used to reduce the (variable-length) lists of nodes and edges into a fixed-length structure that still preserves some properties of the graph. Different instance encodings may have significant impact on performance.","title":"Describing problem instances"},{"location":"usage/#obtaining-heuristic-solutions","text":"By default, LearningSolver uses Machine Learning to accelerate the MIP solution process, while maintaining all optimality guarantees provided by the MIP solver. In the default mode of operation, for example, predicted optimal solutions are used only as MIP starts. For more significant performance benefits, LearningSolver can also be configured to place additional trust in the Machine Learning predictors, by using the mode=\"heuristic\" constructor argument. When operating in this mode, if a ML model is statistically shown (through stratified k-fold cross validation ) to have exceptionally high accuracy, the solver may decide to restrict the search space based on its predictions. The parts of the solution which the ML models cannot predict accurately will still be explored using traditional (branch-and-bound) methods. For particular applications, this mode has been shown to quickly produce optimal or near-optimal solutions (see references and benchmark results ). Danger The heuristic mode provides no optimality guarantees, and therefore should only be used if the solver is first trained on a large and representative set of training instances. Training on a small or non-representative set of instances may produce low-quality solutions, or make the solver incorrectly classify new instances as infeasible.","title":"Obtaining heuristic solutions"},{"location":"usage/#saving-and-loading-solver-state","text":"After solving a large number of training instances, it may be desirable to save the current state of LearningSolver to disk, so that the solver can still use the acquired knowledge after the application restarts. This can be accomplished by using the standard pickle module, as the following example illustrates: from miplearn import LearningSolver import pickle # Solve training instances training_instances = [...] solver = LearningSolver() for instance in training_instances: solver.solve(instance) # Train machine-learning models solver.fit(training_instances) # Save trained solver to disk pickle.dump(solver, open(\"solver.pickle\", \"wb\")) # Application restarts... # Load trained solver from disk solver = pickle.load(open(\"solver.pickle\", \"rb\")) # Solve additional instances test_instances = [...] for instance in test_instances: solver.solve(instance)","title":"Saving and loading solver state"},{"location":"usage/#solving-training-instances-in-parallel","text":"In many situations, training and test instances can be solved in parallel to accelerate the training process. LearningSolver provides the method parallel_solve(instances) to easily achieve this: from miplearn import LearningSolver training_instances = [...] solver = LearningSolver() solver.parallel_solve(training_instances, n_jobs=4) solver.fit(training_instances) # Test phase... test_instances = [...] solver.parallel_solve(test_instances)","title":"Solving training instances in parallel"},{"location":"usage/#current-limitations","text":"Only binary and continuous decision variables are currently supported.","title":"Current Limitations"}]}
\ No newline at end of file
+{"config":{"lang":["en"],"prebuild_index":false,"separator":"[\\s\\-]+"},"docs":[{"location":"","text":"MIPLearn MIPLearn is an extensible framework for Learning-Enhanced Mixed-Integer Optimization , an approach targeted at discrete optimization problems that need to be repeatedly solved with only minor changes to input data. The package uses Machine Learning (ML) to automatically identify patterns in previously solved instances of the problem, or in the solution process itself, and produces hints that can guide a conventional MIP solver towards the optimal solution faster. For particular classes of problems, this approach has been shown to provide significant performance benefits (see benchmark results and references for more details). Features MIPLearn proposes a flexible problem specification format, which allows users to describe their particular optimization problems to a Learning-Enhanced MIP solver, both from the MIP perspective and from the ML perspective, without making any assumptions on the problem being modeled, the mathematical formulation of the problem, or ML encoding. While the format is very flexible, some constraints are enforced to ensure that it is usable by an actual solver. MIPLearn provides a reference implementation of a Learning-Enhanced Solver , which can use the above problem specification format to automatically predict, based on previously solved instances, a number of hints to accelerate MIP performance. Currently, the reference solver is able to predict: (i) partial solutions which are likely to work well as MIP starts; (ii) an initial set of lazy constraints to enforce; (iii) variable branching priorities to accelerate the exploration of the branch-and-bound tree; (iv) the optimal objective value based on the solution to the LP relaxation. The usage of the solver is very straightforward. The most suitable ML models are automatically selected, trained, cross-validated and applied to the problem with no user intervention. MIPLearn provides a set of benchmark problems and random instance generators, covering applications from different domains, which can be used to quickly evaluate new learning-enhanced MIP techniques in a measurable and reproducible way. MIPLearn is customizable and extensible . For MIP and ML researchers exploring new techniques to accelerate MIP performance based on historical data, each component of the reference solver can be individually replaced, extended or customized. Documentation Installation and typical usage Benchmark utilities Benchmark problems, challenges and results Customizing the solver License, authors, references and acknowledgements Souce Code https://github.com/ANL-CEEESA/MIPLearn","title":"Home"},{"location":"#miplearn","text":"MIPLearn is an extensible framework for Learning-Enhanced Mixed-Integer Optimization , an approach targeted at discrete optimization problems that need to be repeatedly solved with only minor changes to input data. The package uses Machine Learning (ML) to automatically identify patterns in previously solved instances of the problem, or in the solution process itself, and produces hints that can guide a conventional MIP solver towards the optimal solution faster. For particular classes of problems, this approach has been shown to provide significant performance benefits (see benchmark results and references for more details).","title":"MIPLearn"},{"location":"#features","text":"MIPLearn proposes a flexible problem specification format, which allows users to describe their particular optimization problems to a Learning-Enhanced MIP solver, both from the MIP perspective and from the ML perspective, without making any assumptions on the problem being modeled, the mathematical formulation of the problem, or ML encoding. While the format is very flexible, some constraints are enforced to ensure that it is usable by an actual solver. MIPLearn provides a reference implementation of a Learning-Enhanced Solver , which can use the above problem specification format to automatically predict, based on previously solved instances, a number of hints to accelerate MIP performance. Currently, the reference solver is able to predict: (i) partial solutions which are likely to work well as MIP starts; (ii) an initial set of lazy constraints to enforce; (iii) variable branching priorities to accelerate the exploration of the branch-and-bound tree; (iv) the optimal objective value based on the solution to the LP relaxation. The usage of the solver is very straightforward. The most suitable ML models are automatically selected, trained, cross-validated and applied to the problem with no user intervention. MIPLearn provides a set of benchmark problems and random instance generators, covering applications from different domains, which can be used to quickly evaluate new learning-enhanced MIP techniques in a measurable and reproducible way. MIPLearn is customizable and extensible . For MIP and ML researchers exploring new techniques to accelerate MIP performance based on historical data, each component of the reference solver can be individually replaced, extended or customized.","title":"Features"},{"location":"#documentation","text":"Installation and typical usage Benchmark utilities Benchmark problems, challenges and results Customizing the solver License, authors, references and acknowledgements","title":"Documentation"},{"location":"#souce-code","text":"https://github.com/ANL-CEEESA/MIPLearn","title":"Souce Code"},{"location":"about/","text":"About Authors Alinson S. Xavier, Argonne National Laboratory < axavier@anl.gov > Feng Qiu, Argonne National Laboratory < fqiu@anl.gov > Acknowledgements Based upon work supported by Laboratory Directed Research and Development (LDRD) funding from Argonne National Laboratory, provided by the Director, Office of Science, of the U.S. Department of Energy under Contract No. DE-AC02-06CH11357. References Learning to Solve Large-Scale Security-Constrained Unit Commitment Problems. Alinson S. Xavier, Feng Qiu, Shabbir Ahmed . INFORMS Journal on Computing (to appear). ArXiv:1902:01696 License MIPLearn, an extensible framework for Learning-Enhanced Mixed-Integer Optimization Copyright \u00a9 2020, UChicago Argonne, LLC. All Rights Reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \"AS IS\" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.","title":"About"},{"location":"about/#about","text":"","title":"About"},{"location":"about/#authors","text":"Alinson S. Xavier, Argonne National Laboratory < axavier@anl.gov > Feng Qiu, Argonne National Laboratory < fqiu@anl.gov >","title":"Authors"},{"location":"about/#acknowledgements","text":"Based upon work supported by Laboratory Directed Research and Development (LDRD) funding from Argonne National Laboratory, provided by the Director, Office of Science, of the U.S. Department of Energy under Contract No. DE-AC02-06CH11357.","title":"Acknowledgements"},{"location":"about/#references","text":"Learning to Solve Large-Scale Security-Constrained Unit Commitment Problems. Alinson S. Xavier, Feng Qiu, Shabbir Ahmed . INFORMS Journal on Computing (to appear). ArXiv:1902:01696","title":"References"},{"location":"about/#license","text":"MIPLearn, an extensible framework for Learning-Enhanced Mixed-Integer Optimization Copyright \u00a9 2020, UChicago Argonne, LLC. All Rights Reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \"AS IS\" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.","title":"License"},{"location":"benchmark/","text":"Benchmarks Utilities Using BenchmarkRunner MIPLearn provides the utility class BenchmarkRunner , which simplifies the task of comparing the performance of different solvers. The snippet below shows its basic usage: from miplearn import BenchmarkRunner, LearningSolver # Create train and test instances train_instances = [...] test_instances = [...] # Training phase... training_solver = LearningSolver(...) training_solver.parallel_solve(train_instances, n_jobs=10) # Test phase... test_solvers = { \"Baseline\": LearningSolver(...), # each solver may have different parameters \"Strategy A\": LearningSolver(...), \"Strategy B\": LearningSolver(...), \"Strategy C\": LearningSolver(...), } benchmark = BenchmarkRunner(test_solvers) benchmark.fit(train_instances) benchmark.parallel_solve(test_instances, n_jobs=2) print(benchmark.raw_results()) The method fit trains the ML models for each individual solver. The method parallel_solve solves the test instances in parallel, and collects solver statistics such as running time and optimal value. Finally, raw_results produces a table of results (Pandas DataFrame) with the following columns: Solver, the name of the solver. Instance, the sequence number identifying the instance. Wallclock Time, the wallclock running time (in seconds) spent by the solver; Lower Bound, the best lower bound obtained by the solver; Upper Bound, the best upper bound obtained by the solver; Gap, the relative MIP integrality gap at the end of the optimization; Nodes, the number of explored branch-and-bound nodes. In addition to the above, there is also a \"Relative\" version of most columns, where the raw number is compared to the solver which provided the best performance. The Relative Wallclock Time for example, indicates how many times slower this run was when compared to the best time achieved by any solver when processing this instance. For example, if this run took 10 seconds, but the fastest solver took only 5 seconds to solve the same instance, the relative wallclock time would be 2. Saving and loading benchmark results When iteratively exploring new formulations, encoding and solver parameters, it is often desirable to avoid repeating parts of the benchmark suite. For example, if the baseline solver has not been changed, there is no need to evaluate its performance again and again when making small changes to the remaining solvers. BenchmarkRunner provides the methods save_results and load_results , which can be used to avoid this repetition, as the next example shows: # Benchmark baseline solvers and save results to a file. benchmark = BenchmarkRunner(baseline_solvers) benchmark.parallel_solve(test_instances) benchmark.save_results(\"baseline_results.csv\") # Benchmark remaining solvers, loading baseline results from file. benchmark = BenchmarkRunner(alternative_solvers) benchmark.load_results(\"baseline_results.csv\") benchmark.fit(training_instances) benchmark.parallel_solve(test_instances)","title":"Benchmark"},{"location":"benchmark/#benchmarks-utilities","text":"","title":"Benchmarks Utilities"},{"location":"benchmark/#using-benchmarkrunner","text":"MIPLearn provides the utility class BenchmarkRunner , which simplifies the task of comparing the performance of different solvers. The snippet below shows its basic usage: from miplearn import BenchmarkRunner, LearningSolver # Create train and test instances train_instances = [...] test_instances = [...] # Training phase... training_solver = LearningSolver(...) training_solver.parallel_solve(train_instances, n_jobs=10) # Test phase... test_solvers = { \"Baseline\": LearningSolver(...), # each solver may have different parameters \"Strategy A\": LearningSolver(...), \"Strategy B\": LearningSolver(...), \"Strategy C\": LearningSolver(...), } benchmark = BenchmarkRunner(test_solvers) benchmark.fit(train_instances) benchmark.parallel_solve(test_instances, n_jobs=2) print(benchmark.raw_results()) The method fit trains the ML models for each individual solver. The method parallel_solve solves the test instances in parallel, and collects solver statistics such as running time and optimal value. Finally, raw_results produces a table of results (Pandas DataFrame) with the following columns: Solver, the name of the solver. Instance, the sequence number identifying the instance. Wallclock Time, the wallclock running time (in seconds) spent by the solver; Lower Bound, the best lower bound obtained by the solver; Upper Bound, the best upper bound obtained by the solver; Gap, the relative MIP integrality gap at the end of the optimization; Nodes, the number of explored branch-and-bound nodes. In addition to the above, there is also a \"Relative\" version of most columns, where the raw number is compared to the solver which provided the best performance. The Relative Wallclock Time for example, indicates how many times slower this run was when compared to the best time achieved by any solver when processing this instance. For example, if this run took 10 seconds, but the fastest solver took only 5 seconds to solve the same instance, the relative wallclock time would be 2.","title":"Using BenchmarkRunner"},{"location":"benchmark/#saving-and-loading-benchmark-results","text":"When iteratively exploring new formulations, encoding and solver parameters, it is often desirable to avoid repeating parts of the benchmark suite. For example, if the baseline solver has not been changed, there is no need to evaluate its performance again and again when making small changes to the remaining solvers. BenchmarkRunner provides the methods save_results and load_results , which can be used to avoid this repetition, as the next example shows: # Benchmark baseline solvers and save results to a file. benchmark = BenchmarkRunner(baseline_solvers) benchmark.parallel_solve(test_instances) benchmark.save_results(\"baseline_results.csv\") # Benchmark remaining solvers, loading baseline results from file. benchmark = BenchmarkRunner(alternative_solvers) benchmark.load_results(\"baseline_results.csv\") benchmark.fit(training_instances) benchmark.parallel_solve(test_instances)","title":"Saving and loading benchmark results"},{"location":"customization/","text":"Customization Selecting the internal MIP solver By default, LearningSolver uses Gurobi as its internal MIP solver. Another supported solver is IBM ILOG CPLEX . To switch between solvers, use the solver constructor argument, as shown below. It is also possible to specify a time limit (in seconds) and a relative MIP gap tolerance. from miplearn import LearningSolver solver = LearningSolver(solver=\"cplex\", time_limit=300, gap_tolerance=1e-3) Selecting solver components LearningSolver is composed by a number of individual machine-learning components, each targeting a different part of the solution process. Each component can be individually enabled, disabled or customized. The following components are enabled by default: LazyConstraintComponent : Predicts which lazy constraint to initially enforce. ObjectiveValueComponent : Predicts the optimal value of the optimization problem, given the optimal solution to the LP relaxation. PrimalSolutionComponent : Predicts optimal values for binary decision variables. In heuristic mode, this component fixes the variables to their predicted values. In exact mode, the predicted values are provided to the solver as a (partial) MIP start. The following components are also available, but not enabled by default: BranchPriorityComponent : Predicts good branch priorities for decision variables. To create a LearningSolver with a specific set of components, the components constructor argument may be used, as the next example shows: # Create a solver without any components solver1 = LearningSolver(components=[]) # Create a solver with only two components solver2 = LearningSolver(components=[ LazyConstraintComponent(...), PrimalSolutionComponent(...), ]) It is also possible to add components to an existing solver using the solver.add method, as shown below. If the solver already holds another component of that type, the new component will replace the previous one. # Create solver with default components solver = LearningSolver() # Replace the default LazyConstraintComponent by one with custom parameters solver.add(LazyConstraintComponent(...)) Adjusting component aggresiveness The aggressiveness of classification components (such as PrimalSolutionComponent and LazyConstraintComponent ) can be adjusted through the threshold constructor argument. Internally, these components ask the ML models how confident they are on each prediction (through the predict_proba method in the sklearn API), and only take into account predictions which have probabilities above the threshold. Lowering a component's threshold increases its aggresiveness, while raising a component's threshold makes it more conservative. MIPLearn also includes MinPrecisionThreshold , a dynamic threshold which adjusts itself automatically during training to achieve a minimum desired true positive rate (also known as precision). The example below shows how to initialize a PrimalSolutionComponent which achieves 95% precision, possibly at the cost of a lower recall. To make the component more aggressive, this precision may be lowered. comp = PrimalSolutionComponent(threshold=MinPrecisionThreshold(0.98))","title":"Customization"},{"location":"customization/#customization","text":"","title":"Customization"},{"location":"customization/#selecting-the-internal-mip-solver","text":"By default, LearningSolver uses Gurobi as its internal MIP solver. Another supported solver is IBM ILOG CPLEX . To switch between solvers, use the solver constructor argument, as shown below. It is also possible to specify a time limit (in seconds) and a relative MIP gap tolerance. from miplearn import LearningSolver solver = LearningSolver(solver=\"cplex\", time_limit=300, gap_tolerance=1e-3)","title":"Selecting the internal MIP solver"},{"location":"customization/#selecting-solver-components","text":"LearningSolver is composed by a number of individual machine-learning components, each targeting a different part of the solution process. Each component can be individually enabled, disabled or customized. The following components are enabled by default: LazyConstraintComponent : Predicts which lazy constraint to initially enforce. ObjectiveValueComponent : Predicts the optimal value of the optimization problem, given the optimal solution to the LP relaxation. PrimalSolutionComponent : Predicts optimal values for binary decision variables. In heuristic mode, this component fixes the variables to their predicted values. In exact mode, the predicted values are provided to the solver as a (partial) MIP start. The following components are also available, but not enabled by default: BranchPriorityComponent : Predicts good branch priorities for decision variables. To create a LearningSolver with a specific set of components, the components constructor argument may be used, as the next example shows: # Create a solver without any components solver1 = LearningSolver(components=[]) # Create a solver with only two components solver2 = LearningSolver(components=[ LazyConstraintComponent(...), PrimalSolutionComponent(...), ]) It is also possible to add components to an existing solver using the solver.add method, as shown below. If the solver already holds another component of that type, the new component will replace the previous one. # Create solver with default components solver = LearningSolver() # Replace the default LazyConstraintComponent by one with custom parameters solver.add(LazyConstraintComponent(...))","title":"Selecting solver components"},{"location":"customization/#adjusting-component-aggresiveness","text":"The aggressiveness of classification components (such as PrimalSolutionComponent and LazyConstraintComponent ) can be adjusted through the threshold constructor argument. Internally, these components ask the ML models how confident they are on each prediction (through the predict_proba method in the sklearn API), and only take into account predictions which have probabilities above the threshold. Lowering a component's threshold increases its aggresiveness, while raising a component's threshold makes it more conservative. MIPLearn also includes MinPrecisionThreshold , a dynamic threshold which adjusts itself automatically during training to achieve a minimum desired true positive rate (also known as precision). The example below shows how to initialize a PrimalSolutionComponent which achieves 95% precision, possibly at the cost of a lower recall. To make the component more aggressive, this precision may be lowered. comp = PrimalSolutionComponent(threshold=MinPrecisionThreshold(0.98))","title":"Adjusting component aggresiveness"},{"location":"problems/","text":"Benchmark Problems, Challenges and Results MIPLearn provides a selection of benchmark problems and random instance generators, covering applications from different fields, that can be used to evaluate new learning-enhanced MIP techniques in a measurable and reproducible way. In this page, we describe these problems, the included instance generators, and we present some benchmark results for LearningSolver with default parameters. Preliminaries Benchmark challenges When evaluating the performance of a conventional MIP solver, benchmark sets , such as MIPLIB and TSPLIB, are typically used. The performance of newly proposed solvers or solution techniques are typically measured as the average (or total) running time the solver takes to solve the entire benchmark set. For Learning-Enhanced MIP solvers, it is also necessary to specify what instances should the solver be trained on (the training instances ) before solving the actual set of instances we are interested in (the test instances ). If the training instances are very similar to the test instances, we would expect a Learning-Enhanced Solver to present stronger perfomance benefits. In MIPLearn, each optimization problem comes with a set of benchmark challenges , which specify how should the training and test instances be generated. The first challenges are typically easier, in the sense that training and test instances are very similar. Later challenges gradually make the sets more distinct, and therefore harder to learn from. Baseline results To illustrate the performance of LearningSolver , and to set a baseline for newly proposed techniques, we present in this page, for each benchmark challenge, a small set of computational results measuring the solution speed of the solver and the solution quality with default parameters. For more detailed computational studies, see references . We compare three solvers: baseline: Gurobi 9.0 with default settings (a conventional state-of-the-art MIP solver) ml-exact: LearningSolver with default settings, using Gurobi 9.0 as internal MIP solver ml-heuristic: Same as above, but with mode=\"heuristic\" All experiments presented here were performed on a Linux server (Ubuntu Linux 18.04 LTS) with Intel Xeon Gold 6230s (2 processors, 40 cores, 80 threads) and 256 GB RAM (DDR4, 2933 MHz). All solvers were restricted to use 4 threads, with no time limits, and 10 instances were solved simultaneously at a time. Maximum Weight Stable Set Problem Problem definition Given a simple undirected graph $G=(V,E)$ and weights $w \\in \\mathbb{R}^V$, the problem is to find a stable set $S \\subseteq V$ that maximizes $ \\sum_{v \\in V} w_v$. We recall that a subset $S \\subseteq V$ is a stable set if no two vertices of $S$ are adjacent. This is one of Karp's 21 NP-complete problems. Random instance generator The class MaxWeightStableSetGenerator can be used to generate random instances of this problem, with user-specified probability distributions. When the constructor parameter fix_graph=True is provided, one random Erd\u0151s-R\u00e9nyi graph $G_{n,p}$ is generated during the constructor, where $n$ and $p$ are sampled from user-provided probability distributions n and p . To generate each instance, the generator independently samples each $w_v$ from the user-provided probability distribution w . When fix_graph=False , a new random graph is generated for each instance, while the remaining parameters are sampled in the same way. Challenge A Fixed random Erd\u0151s-R\u00e9nyi graph $G_{n,p}$ with $n=200$ and $p=5\\%$ Random vertex weights $w_v \\sim U(100, 150)$ 500 training instances, 50 test instances MaxWeightStableSetGenerator(w=uniform(loc=100., scale=50.), n=randint(low=200, high=201), p=uniform(loc=0.05, scale=0.0), fix_graph=True) Traveling Salesman Problem Problem definition Given a list of cities and the distance between each pair of cities, the problem asks for the shortest route starting at the first city, visiting each other city exactly once, then returning to the first city. This problem is a generalization of the Hamiltonian path problem, one of Karp's 21 NP-complete problems. Random problem generator The class TravelingSalesmanGenerator can be used to generate random instances of this problem. Initially, the generator creates $n$ cities $(x_1,y_1),\\ldots,(x_n,y_n) \\in \\mathbb{R}^2$, where $n, x_i$ and $y_i$ are sampled independently from the provided probability distributions n , x and y . For each pair of cities $(i,j)$, the distance $d_{i,j}$ between them is set to: d_{i,j} = \\gamma_{i,j} \\sqrt{(x_i-x_j)^2 + (y_i - y_j)^2} where $\\gamma_{i,j}$ is sampled from the distribution gamma . If fix_cities=True is provided, 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. Challenge A Fixed list of 350 cities in the $[0, 1000]^2$ square $\\gamma_{i,j} \\sim U(0.95, 1.05)$ 500 training instances, 50 test instances 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, round=True, ) Multidimensional 0-1 Knapsack Problem Problem definition Given a set of $n$ items and $m$ types of resources (also called knapsacks ), the problem is to find a subset of items that maximizes profit without consuming more resources than it is available. More precisely, the problem is: \\begin{align*} \\text{maximize} & \\sum_{j=1}^n p_j x_j \\\\ \\text{subject to} & \\sum_{j=1}^n w_{ij} x_j \\leq b_i & \\forall i=1,\\ldots,m \\\\ & x_j \\in \\{0,1\\} & \\forall j=1,\\ldots,n \\end{align*} Random instance generator The class MultiKnapsackGenerator can be used to generate random instances of this problem. The number of items $n$ and knapsacks $m$ are sampled from the user-provided probability distributions n and m . The weights $w_{ij}$ are sampled independently from the provided distribution w . The capacity of knapsack $i$ is set to b_i = \\alpha_i \\sum_{j=1}^n w_{ij} where $\\alpha_i$, the tightness ratio, is sampled from the provided probability distribution alpha . To make the instances more challenging, the costs of the items are linearly correlated to their average weights. More specifically, the price of each item $j$ is set to: p_j = \\sum_{i=1}^m \\frac{w_{ij}}{m} + K u_j, where $K$, the correlation coefficient, and $u_j$, the correlation multiplier, are sampled from the provided probability distributions K and u . If fix_w=True is provided, then $w_{ij}$ are kept the same in all generated instances. This also implies that $n$ and $m$ are kept fixed. Although the prices and capacities are derived from $w_{ij}$, as long as u and K are not constants, the generated instances will still not be completely identical. If a probability distribution w_jitter is provided, then item weights will be set to $w_{ij} \\gamma_{ij}$ where $\\gamma_{ij}$ is sampled from w_jitter . When combined with fix_w=True , this argument may be used to generate instances where the weight of each item is roughly the same, but not exactly identical, across all instances. The prices of the items and the capacities of the knapsacks will be calculated as above, but using these perturbed weights instead. By default, all generated prices, weights and capacities are rounded to the nearest integer number. If round=False is provided, this rounding will be disabled. References Freville, Arnaud, and G\u00e9rard Plateau. An efficient preprocessing procedure for the multidimensional 0\u20131 knapsack problem. Discrete applied mathematics 49.1-3 (1994): 189-212. Fr\u00e9ville, Arnaud. The multidimensional 0\u20131 knapsack problem: An overview. European Journal of Operational Research 155.1 (2004): 1-21. Challenge A 250 variables, 10 constraints, fixed weights $w \\sim U(0, 1000), \\gamma \\sim U(0.95, 1.05)$ $K = 500, u \\sim U(0, 1), \\alpha = 0.25$ 500 training instances, 50 test instances MultiKnapsackGenerator(n=randint(low=250, high=251), m=randint(low=10, high=11), w=uniform(loc=0.0, scale=1000.0), K=uniform(loc=500.0, scale=0.0), u=uniform(loc=0.0, scale=1.0), alpha=uniform(loc=0.25, scale=0.0), fix_w=True, w_jitter=uniform(loc=0.95, scale=0.1), )","title":"Problems"},{"location":"problems/#benchmark-problems-challenges-and-results","text":"MIPLearn provides a selection of benchmark problems and random instance generators, covering applications from different fields, that can be used to evaluate new learning-enhanced MIP techniques in a measurable and reproducible way. In this page, we describe these problems, the included instance generators, and we present some benchmark results for LearningSolver with default parameters.","title":"Benchmark Problems, Challenges and Results"},{"location":"problems/#preliminaries","text":"","title":"Preliminaries"},{"location":"problems/#benchmark-challenges","text":"When evaluating the performance of a conventional MIP solver, benchmark sets , such as MIPLIB and TSPLIB, are typically used. The performance of newly proposed solvers or solution techniques are typically measured as the average (or total) running time the solver takes to solve the entire benchmark set. For Learning-Enhanced MIP solvers, it is also necessary to specify what instances should the solver be trained on (the training instances ) before solving the actual set of instances we are interested in (the test instances ). If the training instances are very similar to the test instances, we would expect a Learning-Enhanced Solver to present stronger perfomance benefits. In MIPLearn, each optimization problem comes with a set of benchmark challenges , which specify how should the training and test instances be generated. The first challenges are typically easier, in the sense that training and test instances are very similar. Later challenges gradually make the sets more distinct, and therefore harder to learn from.","title":"Benchmark challenges"},{"location":"problems/#baseline-results","text":"To illustrate the performance of LearningSolver , and to set a baseline for newly proposed techniques, we present in this page, for each benchmark challenge, a small set of computational results measuring the solution speed of the solver and the solution quality with default parameters. For more detailed computational studies, see references . We compare three solvers: baseline: Gurobi 9.0 with default settings (a conventional state-of-the-art MIP solver) ml-exact: LearningSolver with default settings, using Gurobi 9.0 as internal MIP solver ml-heuristic: Same as above, but with mode=\"heuristic\" All experiments presented here were performed on a Linux server (Ubuntu Linux 18.04 LTS) with Intel Xeon Gold 6230s (2 processors, 40 cores, 80 threads) and 256 GB RAM (DDR4, 2933 MHz). All solvers were restricted to use 4 threads, with no time limits, and 10 instances were solved simultaneously at a time.","title":"Baseline results"},{"location":"problems/#maximum-weight-stable-set-problem","text":"","title":"Maximum Weight Stable Set Problem"},{"location":"problems/#problem-definition","text":"Given a simple undirected graph $G=(V,E)$ and weights $w \\in \\mathbb{R}^V$, the problem is to find a stable set $S \\subseteq V$ that maximizes $ \\sum_{v \\in V} w_v$. We recall that a subset $S \\subseteq V$ is a stable set if no two vertices of $S$ are adjacent. This is one of Karp's 21 NP-complete problems.","title":"Problem definition"},{"location":"problems/#random-instance-generator","text":"The class MaxWeightStableSetGenerator can be used to generate random instances of this problem, with user-specified probability distributions. When the constructor parameter fix_graph=True is provided, one random Erd\u0151s-R\u00e9nyi graph $G_{n,p}$ is generated during the constructor, where $n$ and $p$ are sampled from user-provided probability distributions n and p . To generate each instance, the generator independently samples each $w_v$ from the user-provided probability distribution w . When fix_graph=False , a new random graph is generated for each instance, while the remaining parameters are sampled in the same way.","title":"Random instance generator"},{"location":"problems/#challenge-a","text":"Fixed random Erd\u0151s-R\u00e9nyi graph $G_{n,p}$ with $n=200$ and $p=5\\%$ Random vertex weights $w_v \\sim U(100, 150)$ 500 training instances, 50 test instances MaxWeightStableSetGenerator(w=uniform(loc=100., scale=50.), n=randint(low=200, high=201), p=uniform(loc=0.05, scale=0.0), fix_graph=True)","title":"Challenge A"},{"location":"problems/#traveling-salesman-problem","text":"","title":"Traveling Salesman Problem"},{"location":"problems/#problem-definition_1","text":"Given a list of cities and the distance between each pair of cities, the problem asks for the shortest route starting at the first city, visiting each other city exactly once, then returning to the first city. This problem is a generalization of the Hamiltonian path problem, one of Karp's 21 NP-complete problems.","title":"Problem definition"},{"location":"problems/#random-problem-generator","text":"The class TravelingSalesmanGenerator can be used to generate random instances of this problem. Initially, the generator creates $n$ cities $(x_1,y_1),\\ldots,(x_n,y_n) \\in \\mathbb{R}^2$, where $n, x_i$ and $y_i$ are sampled independently from the provided probability distributions n , x and y . For each pair of cities $(i,j)$, the distance $d_{i,j}$ between them is set to: d_{i,j} = \\gamma_{i,j} \\sqrt{(x_i-x_j)^2 + (y_i - y_j)^2} where $\\gamma_{i,j}$ is sampled from the distribution gamma . If fix_cities=True is provided, 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.","title":"Random problem generator"},{"location":"problems/#challenge-a_1","text":"Fixed list of 350 cities in the $[0, 1000]^2$ square $\\gamma_{i,j} \\sim U(0.95, 1.05)$ 500 training instances, 50 test instances 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, round=True, )","title":"Challenge A"},{"location":"problems/#multidimensional-0-1-knapsack-problem","text":"","title":"Multidimensional 0-1 Knapsack Problem"},{"location":"problems/#problem-definition_2","text":"Given a set of $n$ items and $m$ types of resources (also called knapsacks ), the problem is to find a subset of items that maximizes profit without consuming more resources than it is available. More precisely, the problem is: \\begin{align*} \\text{maximize} & \\sum_{j=1}^n p_j x_j \\\\ \\text{subject to} & \\sum_{j=1}^n w_{ij} x_j \\leq b_i & \\forall i=1,\\ldots,m \\\\ & x_j \\in \\{0,1\\} & \\forall j=1,\\ldots,n \\end{align*}","title":"Problem definition"},{"location":"problems/#random-instance-generator_1","text":"The class MultiKnapsackGenerator can be used to generate random instances of this problem. The number of items $n$ and knapsacks $m$ are sampled from the user-provided probability distributions n and m . The weights $w_{ij}$ are sampled independently from the provided distribution w . The capacity of knapsack $i$ is set to b_i = \\alpha_i \\sum_{j=1}^n w_{ij} where $\\alpha_i$, the tightness ratio, is sampled from the provided probability distribution alpha . To make the instances more challenging, the costs of the items are linearly correlated to their average weights. More specifically, the price of each item $j$ is set to: p_j = \\sum_{i=1}^m \\frac{w_{ij}}{m} + K u_j, where $K$, the correlation coefficient, and $u_j$, the correlation multiplier, are sampled from the provided probability distributions K and u . If fix_w=True is provided, then $w_{ij}$ are kept the same in all generated instances. This also implies that $n$ and $m$ are kept fixed. Although the prices and capacities are derived from $w_{ij}$, as long as u and K are not constants, the generated instances will still not be completely identical. If a probability distribution w_jitter is provided, then item weights will be set to $w_{ij} \\gamma_{ij}$ where $\\gamma_{ij}$ is sampled from w_jitter . When combined with fix_w=True , this argument may be used to generate instances where the weight of each item is roughly the same, but not exactly identical, across all instances. The prices of the items and the capacities of the knapsacks will be calculated as above, but using these perturbed weights instead. By default, all generated prices, weights and capacities are rounded to the nearest integer number. If round=False is provided, this rounding will be disabled. References Freville, Arnaud, and G\u00e9rard Plateau. An efficient preprocessing procedure for the multidimensional 0\u20131 knapsack problem. Discrete applied mathematics 49.1-3 (1994): 189-212. Fr\u00e9ville, Arnaud. The multidimensional 0\u20131 knapsack problem: An overview. European Journal of Operational Research 155.1 (2004): 1-21.","title":"Random instance generator"},{"location":"problems/#challenge-a_2","text":"250 variables, 10 constraints, fixed weights $w \\sim U(0, 1000), \\gamma \\sim U(0.95, 1.05)$ $K = 500, u \\sim U(0, 1), \\alpha = 0.25$ 500 training instances, 50 test instances MultiKnapsackGenerator(n=randint(low=250, high=251), m=randint(low=10, high=11), w=uniform(loc=0.0, scale=1000.0), K=uniform(loc=500.0, scale=0.0), u=uniform(loc=0.0, scale=1.0), alpha=uniform(loc=0.25, scale=0.0), fix_w=True, w_jitter=uniform(loc=0.95, scale=0.1), )","title":"Challenge A"},{"location":"usage/","text":"Usage Installation MIPLearn is mainly written in Python, with some components written in Julia. For this reason, both Python 3.6+ and Julia 1.3+ are required. A mixed-integer solver is also required, and its Python bindings must be properly installed. Supported solvers are CPLEX and Gurobi. Optimization problems currently need to be specified in the Pyomo modeling language. A JuMP interface to the package is currently under development. To install MIPLearn, run the following commands: git clone https://github.com/ANL-CEEESA/MIPLearn.git cd MIPLearn make install After installation, the package miplearn should become available to Python. It can be imported as follows: import miplearn Note To install MIPLearn in another Python environment, switch to that enviroment before running make install . To install the package in development mode, run make develop instead. Using LearningSolver The main class provided by this package is LearningSolver , a learning-enhanced MIP solver which uses information from previously solved instances to accelerate the solution of new instances. The following example shows its basic usage: from miplearn import LearningSolver # List of user-provided instances training_instances = [...] test_instances = [...] # Create solver solver = LearningSolver() # Solve all training instances for instance in training_instances: solver.solve(instance) # Learn from training instances solver.fit(training_instances) # Solve all test instances for instance in test_instances: solver.solve(instance) In this example, we have two lists of user-provided instances: training_instances and test_instances . We start by solving all training instances. Since there is no historical information available at this point, the instances will be processed from scratch, with no ML acceleration. After solving each instance, the solver stores within each instance object the optimal solution, the optimal objective value, and other information that can be used to accelerate future solves. After all training instances are solved, we call solver.fit(training_instances) . This instructs the solver to train all its internal machine-learning models based on the solutions of the (solved) trained instances. Subsequent calls to solver.solve(instance) will automatically use the trained Machine Learning models to accelerate the solution process. Describing problem instances Instances to be solved by LearningSolver must derive from the abstract class miplearn.Instance . The following three abstract methods must be implemented: instance.to_model() , which returns a concrete Pyomo model corresponding to the instance; instance.get_instance_features() , which returns a 1-dimensional Numpy array of (numerical) features describing the entire instance; instance.get_variable_features(var_name, index) , which returns a 1-dimensional array of (numerical) features describing a particular decision variable. The first method is used by LearningSolver to construct a concrete Pyomo model, which will be provided to the internal MIP solver. The second and third methods provide an encoding of the instance, which can be used by the ML models to make predictions. In the knapsack problem, for example, an implementation may decide to provide as instance features the average weights, average prices, number of items and the size of the knapsack. The weight and the price of each individual item could be provided as variable features. See src/python/miplearn/problems/knapsack.py for a concrete example. An optional method which can be implemented is instance.get_variable_category(var_name, index) , which returns a category (a string, an integer or any hashable type) for each decision variable. If two variables have the same category, LearningSolver will use the same internal ML model to predict the values of both variables. By default, all variables belong to the \"default\" category, and therefore only one ML model is used for all variables. If the returned category is None , ML predictors will ignore the variable. It is not necessary to have a one-to-one correspondence between features and problem instances. One important (and deliberate) limitation of MIPLearn, however, is that get_instance_features() must always return arrays of same length for all relevant instances of the problem. Similarly, get_variable_features(var_name, index) must also always return arrays of same length for all variables in each category. It is up to the user to decide how to encode variable-length characteristics of the problem into fixed-length vectors. In graph problems, for example, graph embeddings can be used to reduce the (variable-length) lists of nodes and edges into a fixed-length structure that still preserves some properties of the graph. Different instance encodings may have significant impact on performance. Obtaining heuristic solutions By default, LearningSolver uses Machine Learning to accelerate the MIP solution process, while maintaining all optimality guarantees provided by the MIP solver. In the default mode of operation, for example, predicted optimal solutions are used only as MIP starts. For more significant performance benefits, LearningSolver can also be configured to place additional trust in the Machine Learning predictors, by using the mode=\"heuristic\" constructor argument. When operating in this mode, if a ML model is statistically shown (through stratified k-fold cross validation ) to have exceptionally high accuracy, the solver may decide to restrict the search space based on its predictions. The parts of the solution which the ML models cannot predict accurately will still be explored using traditional (branch-and-bound) methods. For particular applications, this mode has been shown to quickly produce optimal or near-optimal solutions (see references and benchmark results ). Danger The heuristic mode provides no optimality guarantees, and therefore should only be used if the solver is first trained on a large and representative set of training instances. Training on a small or non-representative set of instances may produce low-quality solutions, or make the solver incorrectly classify new instances as infeasible. Saving and loading solver state After solving a large number of training instances, it may be desirable to save the current state of LearningSolver to disk, so that the solver can still use the acquired knowledge after the application restarts. This can be accomplished by using the standard pickle module, as the following example illustrates: from miplearn import LearningSolver import pickle # Solve training instances training_instances = [...] solver = LearningSolver() for instance in training_instances: solver.solve(instance) # Train machine-learning models solver.fit(training_instances) # Save trained solver to disk pickle.dump(solver, open(\"solver.pickle\", \"wb\")) # Application restarts... # Load trained solver from disk solver = pickle.load(open(\"solver.pickle\", \"rb\")) # Solve additional instances test_instances = [...] for instance in test_instances: solver.solve(instance) Solving training instances in parallel In many situations, training and test instances can be solved in parallel to accelerate the training process. LearningSolver provides the method parallel_solve(instances) to easily achieve this: from miplearn import LearningSolver training_instances = [...] solver = LearningSolver() solver.parallel_solve(training_instances, n_jobs=4) solver.fit(training_instances) # Test phase... test_instances = [...] solver.parallel_solve(test_instances) Current Limitations Only binary and continuous decision variables are currently supported.","title":"Usage"},{"location":"usage/#usage","text":"","title":"Usage"},{"location":"usage/#installation","text":"MIPLearn is mainly written in Python, with some components written in Julia. For this reason, both Python 3.6+ and Julia 1.3+ are required. A mixed-integer solver is also required, and its Python bindings must be properly installed. Supported solvers are CPLEX and Gurobi. Optimization problems currently need to be specified in the Pyomo modeling language. A JuMP interface to the package is currently under development. To install MIPLearn, run the following commands: git clone https://github.com/ANL-CEEESA/MIPLearn.git cd MIPLearn make install After installation, the package miplearn should become available to Python. It can be imported as follows: import miplearn Note To install MIPLearn in another Python environment, switch to that enviroment before running make install . To install the package in development mode, run make develop instead.","title":"Installation"},{"location":"usage/#using-learningsolver","text":"The main class provided by this package is LearningSolver , a learning-enhanced MIP solver which uses information from previously solved instances to accelerate the solution of new instances. The following example shows its basic usage: from miplearn import LearningSolver # List of user-provided instances training_instances = [...] test_instances = [...] # Create solver solver = LearningSolver() # Solve all training instances for instance in training_instances: solver.solve(instance) # Learn from training instances solver.fit(training_instances) # Solve all test instances for instance in test_instances: solver.solve(instance) In this example, we have two lists of user-provided instances: training_instances and test_instances . We start by solving all training instances. Since there is no historical information available at this point, the instances will be processed from scratch, with no ML acceleration. After solving each instance, the solver stores within each instance object the optimal solution, the optimal objective value, and other information that can be used to accelerate future solves. After all training instances are solved, we call solver.fit(training_instances) . This instructs the solver to train all its internal machine-learning models based on the solutions of the (solved) trained instances. Subsequent calls to solver.solve(instance) will automatically use the trained Machine Learning models to accelerate the solution process.","title":"Using LearningSolver"},{"location":"usage/#describing-problem-instances","text":"Instances to be solved by LearningSolver must derive from the abstract class miplearn.Instance . The following three abstract methods must be implemented: instance.to_model() , which returns a concrete Pyomo model corresponding to the instance; instance.get_instance_features() , which returns a 1-dimensional Numpy array of (numerical) features describing the entire instance; instance.get_variable_features(var_name, index) , which returns a 1-dimensional array of (numerical) features describing a particular decision variable. The first method is used by LearningSolver to construct a concrete Pyomo model, which will be provided to the internal MIP solver. The second and third methods provide an encoding of the instance, which can be used by the ML models to make predictions. In the knapsack problem, for example, an implementation may decide to provide as instance features the average weights, average prices, number of items and the size of the knapsack. The weight and the price of each individual item could be provided as variable features. See src/python/miplearn/problems/knapsack.py for a concrete example. An optional method which can be implemented is instance.get_variable_category(var_name, index) , which returns a category (a string, an integer or any hashable type) for each decision variable. If two variables have the same category, LearningSolver will use the same internal ML model to predict the values of both variables. By default, all variables belong to the \"default\" category, and therefore only one ML model is used for all variables. If the returned category is None , ML predictors will ignore the variable. It is not necessary to have a one-to-one correspondence between features and problem instances. One important (and deliberate) limitation of MIPLearn, however, is that get_instance_features() must always return arrays of same length for all relevant instances of the problem. Similarly, get_variable_features(var_name, index) must also always return arrays of same length for all variables in each category. It is up to the user to decide how to encode variable-length characteristics of the problem into fixed-length vectors. In graph problems, for example, graph embeddings can be used to reduce the (variable-length) lists of nodes and edges into a fixed-length structure that still preserves some properties of the graph. Different instance encodings may have significant impact on performance.","title":"Describing problem instances"},{"location":"usage/#obtaining-heuristic-solutions","text":"By default, LearningSolver uses Machine Learning to accelerate the MIP solution process, while maintaining all optimality guarantees provided by the MIP solver. In the default mode of operation, for example, predicted optimal solutions are used only as MIP starts. For more significant performance benefits, LearningSolver can also be configured to place additional trust in the Machine Learning predictors, by using the mode=\"heuristic\" constructor argument. When operating in this mode, if a ML model is statistically shown (through stratified k-fold cross validation ) to have exceptionally high accuracy, the solver may decide to restrict the search space based on its predictions. The parts of the solution which the ML models cannot predict accurately will still be explored using traditional (branch-and-bound) methods. For particular applications, this mode has been shown to quickly produce optimal or near-optimal solutions (see references and benchmark results ). Danger The heuristic mode provides no optimality guarantees, and therefore should only be used if the solver is first trained on a large and representative set of training instances. Training on a small or non-representative set of instances may produce low-quality solutions, or make the solver incorrectly classify new instances as infeasible.","title":"Obtaining heuristic solutions"},{"location":"usage/#saving-and-loading-solver-state","text":"After solving a large number of training instances, it may be desirable to save the current state of LearningSolver to disk, so that the solver can still use the acquired knowledge after the application restarts. This can be accomplished by using the standard pickle module, as the following example illustrates: from miplearn import LearningSolver import pickle # Solve training instances training_instances = [...] solver = LearningSolver() for instance in training_instances: solver.solve(instance) # Train machine-learning models solver.fit(training_instances) # Save trained solver to disk pickle.dump(solver, open(\"solver.pickle\", \"wb\")) # Application restarts... # Load trained solver from disk solver = pickle.load(open(\"solver.pickle\", \"rb\")) # Solve additional instances test_instances = [...] for instance in test_instances: solver.solve(instance)","title":"Saving and loading solver state"},{"location":"usage/#solving-training-instances-in-parallel","text":"In many situations, training and test instances can be solved in parallel to accelerate the training process. LearningSolver provides the method parallel_solve(instances) to easily achieve this: from miplearn import LearningSolver training_instances = [...] solver = LearningSolver() solver.parallel_solve(training_instances, n_jobs=4) solver.fit(training_instances) # Test phase... test_instances = [...] solver.parallel_solve(test_instances)","title":"Solving training instances in parallel"},{"location":"usage/#current-limitations","text":"Only binary and continuous decision variables are currently supported.","title":"Current Limitations"}]}
\ No newline at end of file
diff --git a/docs/sitemap.xml b/docs/sitemap.xml
index 22c2e8b..05c6e74 100644
--- a/docs/sitemap.xml
+++ b/docs/sitemap.xml
@@ -1,27 +1,27 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <urlset xmlns="http://www.sitemaps.org/schemas/sitemap/0.9"><url>
      <loc>None</loc>
-     <lastmod>2020-04-14</lastmod>
+     <lastmod>2020-05-05</lastmod>
      <changefreq>daily</changefreq>
     </url><url>
      <loc>None</loc>
-     <lastmod>2020-04-14</lastmod>
+     <lastmod>2020-05-05</lastmod>
      <changefreq>daily</changefreq>
     </url><url>
      <loc>None</loc>
-     <lastmod>2020-04-14</lastmod>
+     <lastmod>2020-05-05</lastmod>
      <changefreq>daily</changefreq>
     </url><url>
      <loc>None</loc>
-     <lastmod>2020-04-14</lastmod>
+     <lastmod>2020-05-05</lastmod>
      <changefreq>daily</changefreq>
     </url><url>
      <loc>None</loc>
-     <lastmod>2020-04-14</lastmod>
+     <lastmod>2020-05-05</lastmod>
      <changefreq>daily</changefreq>
     </url><url>
      <loc>None</loc>
-     <lastmod>2020-04-14</lastmod>
+     <lastmod>2020-05-05</lastmod>
      <changefreq>daily</changefreq>
     </url>
 </urlset>
\ No newline at end of file
diff --git a/docs/sitemap.xml.gz b/docs/sitemap.xml.gz
index 5f85502..b80eaed 100644
Binary files a/docs/sitemap.xml.gz and b/docs/sitemap.xml.gz differ
diff --git a/src/docs/customization.md b/src/docs/customization.md
index def1e94..f635361 100644
--- a/src/docs/customization.md
+++ b/src/docs/customization.md
@@ -44,3 +44,20 @@ solver = LearningSolver()
 # Replace the default LazyConstraintComponent by one with custom parameters 
 solver.add(LazyConstraintComponent(...))
 ```
+
+## Adjusting component aggresiveness
+
+The aggressiveness of classification components (such as `PrimalSolutionComponent` and `LazyConstraintComponent`) can
+be adjusted through the `threshold` constructor argument. Internally, these components ask the ML models how confident
+they are on each prediction (through the `predict_proba` method in the sklearn API), and only take into account
+predictions which have probabilities above the threshold. Lowering a component's threshold increases its aggresiveness,
+while raising a component's threshold makes it more conservative. 
+
+MIPLearn also includes `MinPrecisionThreshold`, a dynamic threshold which adjusts itself automatically during training
+to achieve a minimum desired true positive rate (also known as precision). The example below shows how to initialize
+a `PrimalSolutionComponent` which achieves 95% precision, possibly at the cost of a lower recall. To make the component
+more aggressive, this precision may be lowered.
+
+```python
+comp = PrimalSolutionComponent(threshold=MinPrecisionThreshold(0.98))
+```
diff --git a/src/python/miplearn/__init__.py b/src/python/miplearn/__init__.py
index 98d5adb..265adb6 100644
--- a/src/python/miplearn/__init__.py
+++ b/src/python/miplearn/__init__.py
@@ -11,7 +11,7 @@ from .components.component import Component
 from .components.objective import ObjectiveValueComponent
 from .components.lazy import LazyConstraintsComponent
 from .components.primal import PrimalSolutionComponent
-from .components.branching import BranchPriorityComponent
+from .components.branching import BranchPriorityComponent, BranchPriorityExtractor
 
 from .classifiers.adaptive import AdaptiveClassifier
 
diff --git a/src/python/miplearn/classifiers/adaptive.py b/src/python/miplearn/classifiers/adaptive.py
index beb41ac..2f6bf8a 100644
--- a/src/python/miplearn/classifiers/adaptive.py
+++ b/src/python/miplearn/classifiers/adaptive.py
@@ -5,106 +5,62 @@
 import logging
 from copy import deepcopy
 
-import numpy as np
-
 from miplearn.classifiers import Classifier
-from sklearn.model_selection import cross_val_score
+from miplearn.classifiers.counting import CountingClassifier
+from miplearn.classifiers.evaluator import ClassifierEvaluator
+from sklearn.linear_model import LogisticRegression
+from sklearn.neighbors import KNeighborsClassifier
+from sklearn.pipeline import make_pipeline
+from sklearn.preprocessing import StandardScaler
 
 logger = logging.getLogger(__name__)
 
 
 class AdaptiveClassifier(Classifier):
     """
-    A classifier that automatically switches strategies based on the number of
-    samples and cross-validation scores.
+    A meta-classifier which dynamically selects what actual classifier to use
+    based on its cross-validation score on a particular training data set.
     """
+
     def __init__(self,
-                 predictor=None,
-                 min_samples_predict=1,
-                 min_samples_cv=100,
-                 thr_fix=0.999,
-                 thr_alpha=0.50,
-                 thr_balance=0.95,
-                 ):
-        self.min_samples_predict = min_samples_predict
-        self.min_samples_cv = min_samples_cv
-        self.thr_fix = thr_fix
-        self.thr_alpha = thr_alpha
-        self.thr_balance = thr_balance
-        self.predictor_factory = predictor
-        self.predictor = None
+                 candidates=None,
+                 evaluator=ClassifierEvaluator()):
+        """
+        Initializes the meta-classifier.
+        """
+        if candidates is None:
+            candidates = {
+                "knn(100)": {
+                    "classifier": KNeighborsClassifier(n_neighbors=100),
+                    "min samples": 100,
+                },
+                "logistic": {
+                    "classifier": make_pipeline(StandardScaler(),
+                                                LogisticRegression()),
+                    "min samples": 30,
+                },
+                "counting": {
+                    "classifier": CountingClassifier(),
+                    "min samples": 0,
+                }
+            }
+        self.candidates = candidates
+        self.evaluator = evaluator
+        self.classifier = None
 
     def fit(self, x_train, y_train):
+        best_name, best_clf, best_score = None, None, -float("inf")
         n_samples = x_train.shape[0]
-
-        # If number of samples is too small, don't predict anything.
-        if n_samples < self.min_samples_predict:
-            logger.debug("    Too few samples (%d); always predicting false" % n_samples)
-            self.predictor = 0
-            return
-
-        # If vast majority of observations are false, always return false.
-        y_train_avg = np.average(y_train)
-        if y_train_avg <= 1.0 - self.thr_fix:
-            logger.debug("    Most samples are negative (%.3f); always returning false" % y_train_avg)
-            self.predictor = 0
-            return
-
-        # If vast majority of observations are true, always return true.
-        if y_train_avg >= self.thr_fix:
-            logger.debug("    Most samples are positive (%.3f); always returning true" % y_train_avg)
-            self.predictor = 1
-            return
-
-        # If classes are too unbalanced, don't predict anything.
-        if y_train_avg < (1 - self.thr_balance) or y_train_avg > self.thr_balance:
-            logger.debug("    Classes are too unbalanced (%.3f); always returning false" % y_train_avg)
-            self.predictor = 0
-            return
-
-        # Select ML model if none is provided
-        if self.predictor_factory is None:
-            if n_samples < 30:
-                from sklearn.neighbors import KNeighborsClassifier
-                self.predictor_factory = KNeighborsClassifier(n_neighbors=n_samples)
-            else:
-                from sklearn.pipeline import make_pipeline
-                from sklearn.preprocessing import StandardScaler
-                from sklearn.linear_model import LogisticRegression
-                self.predictor_factory = make_pipeline(StandardScaler(), LogisticRegression())
-
-        # Create predictor
-        if callable(self.predictor_factory):
-            pred = self.predictor_factory()
-        else:
-            pred = deepcopy(self.predictor_factory)
-
-        # Skip cross-validation if number of samples is too small
-        if n_samples < self.min_samples_cv:
-            logger.debug("    Too few samples (%d); skipping cross validation" % n_samples)
-            self.predictor = pred
-            self.predictor.fit(x_train, y_train)
-            return
-
-        # Calculate cross-validation score
-        cv_score = np.mean(cross_val_score(pred, x_train, y_train, cv=5))
-        dummy_score = max(y_train_avg, 1 - y_train_avg)
-        cv_thr = 1. * self.thr_alpha + dummy_score * (1 - self.thr_alpha)
-
-        # If cross-validation score is too low, don't predict anything.
-        if cv_score < cv_thr:
-            logger.debug("    Score is too low (%.3f < %.3f); always returning false" % (cv_score, cv_thr))
-            self.predictor = 0
-        else:
-            logger.debug("    Score is acceptable (%.3f > %.3f); training classifier" % (cv_score, cv_thr))
-            self.predictor = pred
-            self.predictor.fit(x_train, y_train)
+        for (name, clf_dict) in self.candidates.items():
+            if n_samples < clf_dict["min samples"]:
+                continue
+            clf = deepcopy(clf_dict["classifier"])
+            clf.fit(x_train, y_train)
+            score = self.evaluator.evaluate(clf, x_train, y_train)
+            if score > best_score:
+                best_name, best_clf, best_score = name, clf, score
+        logger.debug("Best classifier: %s (score=%.3f)" % (best_name, best_score))
+        self.classifier = best_clf
 
     def predict_proba(self, x_test):
-        if isinstance(self.predictor, int):
-            y_pred = np.zeros((x_test.shape[0], 2))
-            y_pred[:, self.predictor] = 1.0
-            return y_pred
-        else:
-            return self.predictor.predict_proba(x_test)
-
+        return self.classifier.predict_proba(x_test)
diff --git a/src/python/miplearn/classifiers/counting.py b/src/python/miplearn/classifiers/counting.py
index 4b6508a..51cb0bf 100644
--- a/src/python/miplearn/classifiers/counting.py
+++ b/src/python/miplearn/classifiers/counting.py
@@ -21,7 +21,8 @@ class CountingClassifier(Classifier):
         self.mean = np.mean(y_train)
 
     def predict_proba(self, x_test):
-        return np.array([[1 - self.mean, self.mean]])
+        return np.array([[1 - self.mean, self.mean]
+                         for _ in range(x_test.shape[0])])
 
     def __repr__(self):
-        return "CountingClassifier(mean=%.3f)" % self.mean
+        return "CountingClassifier(mean=%s)" % self.mean
diff --git a/src/python/miplearn/classifiers/evaluator.py b/src/python/miplearn/classifiers/evaluator.py
new file mode 100644
index 0000000..5ff0f30
--- /dev/null
+++ b/src/python/miplearn/classifiers/evaluator.py
@@ -0,0 +1,15 @@
+#  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 sklearn.metrics import roc_auc_score
+
+
+class ClassifierEvaluator:
+    def __init__(self):
+        pass
+
+    def evaluate(self, clf, x_train, y_train):
+        # FIXME: use cross-validation
+        proba = clf.predict_proba(x_train)
+        return roc_auc_score(y_train, proba[:, 1])
diff --git a/src/python/miplearn/classifiers/tests/test_counting.py b/src/python/miplearn/classifiers/tests/test_counting.py
index d7ad627..0f3e6f1 100644
--- a/src/python/miplearn/classifiers/tests/test_counting.py
+++ b/src/python/miplearn/classifiers/tests/test_counting.py
@@ -12,6 +12,7 @@ E = 0.1
 def test_counting():
     clf = CountingClassifier()
     clf.fit(np.zeros((8, 25)), [0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0])
-    expected_proba = np.array([[0.375, 0.625]])
-    actual_proba = clf.predict_proba(np.zeros((1, 25)))
+    expected_proba = np.array([[0.375, 0.625],
+                               [0.375, 0.625]])
+    actual_proba = clf.predict_proba(np.zeros((2, 25)))
     assert norm(actual_proba - expected_proba) < E
diff --git a/src/python/miplearn/classifiers/tests/test_evaluator.py b/src/python/miplearn/classifiers/tests/test_evaluator.py
new file mode 100644
index 0000000..a7f16d6
--- /dev/null
+++ b/src/python/miplearn/classifiers/tests/test_evaluator.py
@@ -0,0 +1,20 @@
+#  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.
+
+import numpy as np
+from miplearn.classifiers.evaluator import ClassifierEvaluator
+from sklearn.neighbors import KNeighborsClassifier
+
+
+def test_evaluator():
+    clf_a = KNeighborsClassifier(n_neighbors=1)
+    clf_b = KNeighborsClassifier(n_neighbors=2)
+    x_train = np.array([[0, 0], [1, 0]])
+    y_train = np.array([0, 1])
+    clf_a.fit(x_train, y_train)
+    clf_b.fit(x_train, y_train)
+    ev = ClassifierEvaluator()
+    assert ev.evaluate(clf_a, x_train, y_train) == 1.0
+    assert ev.evaluate(clf_b, x_train, y_train) == 0.5
+
diff --git a/src/python/miplearn/classifiers/tests/test_threshold.py b/src/python/miplearn/classifiers/tests/test_threshold.py
new file mode 100644
index 0000000..6d224b1
--- /dev/null
+++ b/src/python/miplearn/classifiers/tests/test_threshold.py
@@ -0,0 +1,34 @@
+#  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 unittest.mock import Mock
+
+import numpy as np
+from miplearn.classifiers import Classifier
+from miplearn.classifiers.threshold import MinPrecisionThreshold
+
+
+def test_threshold_dynamic():
+    clf = Mock(spec=Classifier)
+    clf.predict_proba = Mock(return_value=np.array([
+        [0.10, 0.90],
+        [0.10, 0.90],
+        [0.20, 0.80],
+        [0.30, 0.70],
+    ]))
+    x_train = np.array([0, 1, 2, 3])
+    y_train = np.array([1, 1, 0, 0])
+
+    threshold = MinPrecisionThreshold(min_precision=1.0)
+    assert threshold.find(clf, x_train, y_train) == 0.90
+
+    threshold = MinPrecisionThreshold(min_precision=0.65)
+    assert threshold.find(clf, x_train, y_train) == 0.80
+
+    threshold = MinPrecisionThreshold(min_precision=0.50)
+    assert threshold.find(clf, x_train, y_train) == 0.70
+
+    threshold = MinPrecisionThreshold(min_precision=0.00)
+    assert threshold.find(clf, x_train, y_train) == 0.70
+
diff --git a/src/python/miplearn/classifiers/threshold.py b/src/python/miplearn/classifiers/threshold.py
new file mode 100644
index 0000000..2026506
--- /dev/null
+++ b/src/python/miplearn/classifiers/threshold.py
@@ -0,0 +1,45 @@
+#  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 abc import abstractmethod, ABC
+
+import numpy as np
+from sklearn.metrics._ranking import _binary_clf_curve
+
+
+class DynamicThreshold(ABC):
+    @abstractmethod
+    def find(self, clf, x_train, y_train):
+        """
+        Given a trained binary classifier `clf` and a training data set,
+        returns the numerical threshold (float) satisfying some criterea.
+        """
+        pass
+
+
+class MinPrecisionThreshold(DynamicThreshold):
+    """
+    The smallest possible threshold satisfying a minimum acceptable true
+    positive rate (also known as precision).
+    """
+
+    def __init__(self, min_precision):
+        self.min_precision = min_precision
+
+    def find(self, clf, x_train, y_train):
+        proba = clf.predict_proba(x_train)
+
+        assert isinstance(proba, np.ndarray), \
+            "classifier should return numpy array"
+        assert proba.shape == (x_train.shape[0], 2), \
+            "classifier should return (%d,%d)-shaped array, not %s" % (
+                x_train.shape[0], 2, str(proba.shape))
+
+        fps, tps, thresholds = _binary_clf_curve(y_train, proba[:, 1])
+        precision = tps / (tps + fps)
+
+        for k in reversed(range(len(precision))):
+            if precision[k] >= self.min_precision:
+                return thresholds[k]
+        return 2.0
diff --git a/src/python/miplearn/components/branching.py b/src/python/miplearn/components/branching.py
index da66b25..5979868 100644
--- a/src/python/miplearn/components/branching.py
+++ b/src/python/miplearn/components/branching.py
@@ -2,147 +2,119 @@
 #  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 Extractor
-from abc import ABC, abstractmethod
-from sklearn.neighbors import KNeighborsRegressor
+import logging
+import os
+import subprocess
+import tempfile
+from copy import deepcopy
+
 import numpy as np
-from tqdm.auto import tqdm
-from joblib import Parallel, delayed
-import multiprocessing
+from pyomo.core import Var
+from pyomo.core.base.label import TextLabeler
+from sklearn.neighbors import KNeighborsRegressor
+from tqdm import tqdm
 
+from .component import Component
+from ..extractors import Extractor, VariableFeaturesExtractor
 
-def _default_branch_priority_predictor():
-    return KNeighborsRegressor(n_neighbors=1)
+logger = logging.getLogger(__name__)
+
+
+class BranchPriorityExtractor(Extractor):
+    def extract(self, instances):
+        result = {}
+        for instance in tqdm(instances,
+                             desc="Extract (branch)",
+                             disable=len(instances) < 5):
+            var_split = self.split_variables(instance)
+            for (category, var_index_pairs) in var_split.items():
+                if category not in result:
+                    result[category] = []
+                for (var_name, index) in var_index_pairs:
+                    result[category] += [instance.branch_priorities[var_name][index]]
+        for category in result:
+            result[category] = np.array(result[category])
+        return result
 
 
 class BranchPriorityComponent(Component):
     def __init__(self,
                  node_limit=10_000,
-                 predictor=_default_branch_priority_predictor,
-                ):
-        self.pending_instances = []
-        self.x_train = {}
-        self.y_train = {}
-        self.predictors = {}
+                 regressor=KNeighborsRegressor(n_neighbors=1),
+                 ):
         self.node_limit = node_limit
-        self.predictor_factory = predictor
-    
+        self.regressors = {}
+        self.regressor_prototype = regressor
+
     def before_solve(self, solver, instance, model):
-        assert solver.internal_solver.name == "gurobi_persistent", "Only GurobiPersistent is currently supported"
-        from gurobipy import GRB
-        var_split = Extractor.split_variables(instance, model)
-        for category in var_split.keys():
-            if category not in self.predictors.keys():
-                continue
-            var_index_pairs = var_split[category]
-            for (i, (var, index)) in enumerate(var_index_pairs):
-                x = self._build_x(instance, var, index)
-                y = self.predictors[category].predict([x])[0][0]
-                gvar = solver.internal_solver._pyomo_var_to_solver_var_map[var[index]]
-                gvar.setAttr(GRB.Attr.BranchPriority, int(round(y)))
-
-                    
-    def after_solve(self, solver, instance, model):
-        self.pending_instances += [instance]
-    
-    def fit(self, solver, n_jobs=1):
-        def _process(instance):
-            # Create LP file
-            import subprocess, tempfile, os, sys
-            lp_file = tempfile.NamedTemporaryFile(suffix=".lp")
-            priority_file = tempfile.NamedTemporaryFile()
+        logger.info("Predicting branching priorities...")
+        priorities = self.predict(instance)
+        solver.internal_solver.set_branching_priorities(priorities)
+
+    def after_solve(self, solver, instance, model, results):
+        pass
+
+    def fit(self, training_instances, n_jobs=1):
+        for instance in tqdm(training_instances, desc="Fit (branch)"):
+            if not hasattr(instance, "branch_priorities"):
+                instance.branch_priorities = self.compute_priorities(instance)
+        x, y = self.x(training_instances), self.y(training_instances)
+        for category in x.keys():
+            self.regressors[category] = deepcopy(self.regressor_prototype)
+            self.regressors[category].fit(x[category], y[category])
+
+    def x(self, instances):
+        return VariableFeaturesExtractor().extract(instances)
+
+    def y(self, instances):
+        return BranchPriorityExtractor().extract(instances)
+
+    def compute_priorities(self, instance, model=None):
+        # Create LP file
+        lp_file = tempfile.NamedTemporaryFile(suffix=".lp")
+        if model is None:
             model = instance.to_model()
-            model.write(lp_file.name)
-            
-            # Run Julia script
-            src_dirname = os.path.dirname(os.path.realpath(__file__))
-            julia_dirname = "%s/../../../julia" % src_dirname
-            priority_file = tempfile.NamedTemporaryFile(mode="r")
-            subprocess.run(["julia",
-                            "--project=%s" % julia_dirname,
-                            "%s/src/branching.jl" % julia_dirname,
-                            lp_file.name,
-                            priority_file.name,
-                            str(self.node_limit),
-                           ],
-                           check=True,
-                          )
-            
-            # Parse output
-            tokens = [line.strip().split(",") for line in priority_file.readlines()]
-            lp_varname_to_priority = {t[0]: int(t[1]) for t in tokens}
-        
-            # Map priorities back to Pyomo variables
-            pyomo_var_to_priority = {}
-            from pyomo.core import Var
-            from pyomo.core.base.label import TextLabeler
-            labeler = TextLabeler()
-            symbol_map = list(model.solutions.symbol_map.values())[0]
-            
-            # Build x_train and y_train
-            comp = BranchPriorityComponent()
-            for var in model.component_objects(Var):
-                for index in var:
-                    category = instance.get_variable_category(var, index)
-                    if category is None:
-                        continue
-                    lp_varname = symbol_map.getSymbol(var[index], labeler)
-                    var_priority = lp_varname_to_priority[lp_varname]
-                    x = self._build_x(instance, var, index)
-                    y = np.array([var_priority])
-                    
-                    if category not in comp.x_train.keys():
-                        comp.x_train[category] = np.array([x])
-                        comp.y_train[category] = np.array([y])
-                    else:
-                        comp.x_train[category] = np.vstack([comp.x_train[category], x])
-                        comp.y_train[category] = np.vstack([comp.y_train[category], y])
-                
-            return comp
-        
-        # Run strong branching on pending instances
-        subcomponents = Parallel(n_jobs=n_jobs)(
-            delayed(_process)(instance)
-            for instance in tqdm(self.pending_instances, desc="Fit (branch)")
-        )
-        self.merge(subcomponents)
-        self.pending_instances.clear()
-
-        # Retrain ML predictors
-        for category in self.x_train.keys():
-            x_train = self.x_train[category]
-            y_train = self.y_train[category]
-            self.predictors[category] = self.predictor_factory()
-            self.predictors[category].fit(x_train, y_train)
-        
-    def _build_x(self, instance, var, index):
-        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]
\ No newline at end of file
+        model.write(lp_file.name)
+
+        # Run Julia script
+        src_dirname = os.path.dirname(os.path.realpath(__file__))
+        julia_dirname = "%s/../../../julia" % src_dirname
+        priority_file = tempfile.NamedTemporaryFile(mode="r")
+        subprocess.run(["julia",
+                        "--project=%s" % julia_dirname,
+                        "%s/src/branching.jl" % julia_dirname,
+                        lp_file.name,
+                        priority_file.name,
+                        str(self.node_limit)],
+                       check=True)
+
+        # Parse output
+        tokens = [line.strip().split(",") for line in priority_file.readlines()]
+        lp_varname_to_priority = {t[0]: int(t[1]) for t in tokens}
+
+        # Map priorities back to Pyomo variables
+        labeler = TextLabeler()
+        symbol_map = list(model.solutions.symbol_map.values())[0]
+        priorities = {}
+        for var in model.component_objects(Var):
+            priorities[var.name] = {}
+            for index in var:
+                category = instance.get_variable_category(var, index)
+                if category is None:
+                    continue
+                lp_varname = symbol_map.getSymbol(var[index], labeler)
+                var_priority = lp_varname_to_priority[lp_varname]
+                priorities[var.name][index] = var_priority
+        return priorities
+
+    def predict(self, instance):
+        priority = {}
+        x_test = self.x([instance])
+        var_split = Extractor.split_variables(instance)
+        for category in self.regressors.keys():
+            y_test = self.regressors[category].predict(x_test[category])
+            for (i, (var, index)) in enumerate(var_split[category]):
+                if var not in priority.keys():
+                    priority[var] = {}
+                priority[var][index] = y_test[i]
+        return priority
diff --git a/src/python/miplearn/components/primal.py b/src/python/miplearn/components/primal.py
index 1ce3a7f..d430768 100644
--- a/src/python/miplearn/components/primal.py
+++ b/src/python/miplearn/components/primal.py
@@ -4,12 +4,10 @@
 
 from copy import deepcopy
 
-from miplearn.classifiers.adaptive import AdaptiveClassifier
-from miplearn.components import classifier_evaluation_dict
-from sklearn.metrics import roc_curve
-from p_tqdm import p_map
-
 from .component import Component
+from ..classifiers.adaptive import AdaptiveClassifier
+from ..classifiers.threshold import MinPrecisionThreshold, DynamicThreshold
+from ..components import classifier_evaluation_dict
 from ..extractors import *
 
 logger = logging.getLogger(__name__)
@@ -23,18 +21,12 @@ class PrimalSolutionComponent(Component):
     def __init__(self,
                  classifier=AdaptiveClassifier(),
                  mode="exact",
-                 max_fpr=[1e-3, 1e-3],
-                 min_threshold=[0.75, 0.75],
-                 dynamic_thresholds=True,
-                 ):
+                 threshold=MinPrecisionThreshold(0.98)):
         self.mode = mode
-        self.is_warm_start_available = False
-        self.max_fpr = max_fpr
-        self.min_threshold = min_threshold
-        self.thresholds = {}
         self.classifiers = {}
+        self.thresholds = {}
+        self.threshold_prototype = threshold
         self.classifier_prototype = classifier
-        self.dynamic_thresholds = dynamic_thresholds
 
     def before_solve(self, solver, instance, model):
         solution = self.predict(instance)
@@ -46,79 +38,64 @@ class PrimalSolutionComponent(Component):
     def after_solve(self, solver, instance, model, results):
         pass
 
+    def x(self, training_instances):
+        return VariableFeaturesExtractor().extract(training_instances)
+
+    def y(self, training_instances):
+        return SolutionExtractor().extract(training_instances)
+
     def fit(self, training_instances, n_jobs=1):
         logger.debug("Extracting features...")
         features = VariableFeaturesExtractor().extract(training_instances)
         solutions = SolutionExtractor().extract(training_instances)
 
-        def _fit(args):
-            category, label = args[0], args[1]
+        for category in tqdm(features.keys(), desc="Fit (primal)"):
             x_train = features[category]
-            y_train = solutions[category]
-            y = y_train[:, label].astype(int)
-
-            if isinstance(self.classifier_prototype, list):
-                clf = deepcopy(self.classifier_prototype[label])
-            else:
-                clf = deepcopy(self.classifier_prototype)
-            clf.fit(x_train, y)
-
-            y_avg = np.average(y)
-            if (not self.dynamic_thresholds) or y_avg <= 0.001 or y_avg >= 0.999:
-                return {"classifier": clf,
-                        "threshold": self.min_threshold[label]}
-
-            proba = clf.predict_proba(x_train)
-            assert isinstance(proba, np.ndarray), \
-                "classifier should return numpy array"
-            assert proba.shape == (x_train.shape[0], 2), \
-                "classifier should return (%d,%d)-shaped array, not %s" % (
-                    x_train.shape[0], 2, str(proba.shape))
-
-            y_scores = proba[:, 1]
-            fpr, tpr, thresholds = roc_curve(y, y_scores)
-            k = 0
-            while True:
-                if (k + 1) > len(fpr):
-                    break
-                if fpr[k + 1] > self.max_fpr[label]:
-                    break
-                if thresholds[k + 1] < self.min_threshold[label]:
-                    break
-                k = k + 1
-            self.thresholds[category, label] = thresholds[k]
-
-            return {"classifier": clf,
-                    "threshold": thresholds[k]}
-
-        items = [(category, label)
-                 for category in features.keys()
-                 for label in [0, 1]]
-
-        if n_jobs == 1:
-            results = list(map(_fit, tqdm(items, desc="Fit (primal)")))
-        else:
-            results = p_map(_fit, items, num_cpus=n_jobs)
-
-        for (idx, (category, label)) in enumerate(items):
-            self.thresholds[category, label] = results[idx]["threshold"]
-            self.classifiers[category, label] = results[idx]["classifier"]
+            for label in [0, 1]:
+                y_train = solutions[category][:, label].astype(int)
+
+                # If all samples are either positive or negative, make constant predictions
+                y_avg = np.average(y_train)
+                if y_avg < 0.001 or y_avg >= 0.999:
+                    self.classifiers[category, label] = round(y_avg)
+                    self.thresholds[category, label] = 0.50
+                    continue
+
+                # Create a copy of classifier prototype and train it
+                if isinstance(self.classifier_prototype, list):
+                    clf = deepcopy(self.classifier_prototype[label])
+                else:
+                    clf = deepcopy(self.classifier_prototype)
+                clf.fit(x_train, y_train)
+
+                # Find threshold (dynamic or static)
+                if isinstance(self.threshold_prototype, DynamicThreshold):
+                    self.thresholds[category, label] = self.threshold_prototype.find(clf, x_train, y_train)
+                else:
+                    self.thresholds[category, label] = deepcopy(self.threshold_prototype)
+
+                self.classifiers[category, label] = clf
 
     def predict(self, instance):
-        x_test = VariableFeaturesExtractor().extract([instance])
         solution = {}
+        x_test = VariableFeaturesExtractor().extract([instance])
         var_split = Extractor.split_variables(instance)
         for category in var_split.keys():
+            n = len(var_split[category])
             for (i, (var, index)) in enumerate(var_split[category]):
                 if var not in solution.keys():
                     solution[var] = {}
                 solution[var][index] = None
-                for label in [0, 1]:
-                    if (category, label) not in self.classifiers.keys():
-                        continue
-                    ws = self.classifiers[category, label].predict_proba(x_test[category])
-                    logger.debug("%s[%s] ws=%.6f threshold=%.6f" %
-                                 (var, index, ws[i, 1], self.thresholds[category, label]))
+            for label in [0, 1]:
+                if (category, label) not in self.classifiers.keys():
+                    continue
+                clf = self.classifiers[category, label]
+                if isinstance(clf, float):
+                    ws = np.array([[1 - clf, clf] for _ in range(n)])
+                else:
+                    ws = clf.predict_proba(x_test[category])
+                assert ws.shape == (n, 2), "ws.shape should be (%d, 2) not %s" % (n, ws.shape)
+                for (i, (var, index)) in enumerate(var_split[category]):
                     if ws[i, 1] >= self.thresholds[category, label]:
                         solution[var][index] = label
         return solution
diff --git a/src/python/miplearn/components/tests/test_branching.py b/src/python/miplearn/components/tests/test_branching.py
index b17f4f7..c0a264b 100644
--- a/src/python/miplearn/components/tests/test_branching.py
+++ b/src/python/miplearn/components/tests/test_branching.py
@@ -1,49 +1,45 @@
 #  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 unittest.mock import Mock
 
-from miplearn import BranchPriorityComponent, LearningSolver
-from miplearn.problems.knapsack import KnapsackInstance
 import numpy as np
-import tempfile
-
-def _get_instances():
-    return [
-        KnapsackInstance(
-            weights=[23., 26., 20., 18.],
-            prices=[505., 352., 458., 220.],
-            capacity=67.,
-        ),
-    ] * 2
-
-
-# def test_branching():
-#     instances = _get_instances()
-#     component = BranchPriorityComponent()
-#     for instance in instances:
-#         component.after_solve(None, instance, None)
-#     component.fit(None)
-#     for key in ["default"]:
-#         assert key in component.x_train.keys()
-#         assert key in component.y_train.keys()
-#         assert component.x_train[key].shape == (8,  4)
-#         assert component.y_train[key].shape == (8,  1)
-        
-        
-# def test_branch_priority_save_load():
-#     state_file = tempfile.NamedTemporaryFile(mode="r")
-#     solver = LearningSolver(components={"branch-priority": BranchPriorityComponent()})
-#     solver.parallel_solve(_get_instances(), n_jobs=2)
-#     solver.fit()
-#     comp = solver.components["branch-priority"]
-#     assert comp.x_train["default"].shape == (8, 4)
-#     assert comp.y_train["default"].shape == (8, 1)
-#     assert "default" in comp.predictors.keys()
-#     solver.save_state(state_file.name)
-#
-#     solver = LearningSolver(components={"branch-priority": BranchPriorityComponent()})
-#     solver.load_state(state_file.name)
-#     comp = solver.components["branch-priority"]
-#     assert comp.x_train["default"].shape == (8, 4)
-#     assert comp.y_train["default"].shape == (8, 1)
-#     assert "default" in comp.predictors.keys()
+from miplearn import BranchPriorityComponent, BranchPriorityExtractor
+from miplearn.classifiers import Regressor
+from miplearn.tests import get_training_instances_and_models
+
+
+def test_branch_extract():
+    instances, models = get_training_instances_and_models()
+    instances[0].branch_priorities = {"x": {0: 100, 1: 200, 2: 300, 3: 400}}
+    instances[1].branch_priorities = {"x": {0: 150, 1: 250, 2: 350, 3: 450}}
+    priorities = BranchPriorityExtractor().extract(instances)
+    assert priorities["default"].tolist() == [100, 200, 300, 400, 150, 250, 350, 450]
+
+
+def test_branch_calculate():
+    instances, models = get_training_instances_and_models()
+    comp = BranchPriorityComponent()
+
+    # If instances do not have branch_priority property, fit should compute them
+    comp.fit(instances)
+    assert instances[0].branch_priorities == {"x": {0: 5730, 1: 24878, 2: 0, 3: 0,}}
+
+    # If instances already have branch_priority, fit should not modify them
+    instances[0].branch_priorities = {"x": {0: 100, 1: 200, 2: 300, 3: 400}}
+    comp.fit(instances)
+    assert instances[0].branch_priorities == {"x": {0: 100, 1: 200, 2: 300, 3: 400}}
+
+
+def test_branch_x_y_predict():
+    instances, models = get_training_instances_and_models()
+    instances[0].branch_priorities = {"x": {0: 100, 1: 200, 2: 300, 3: 400}}
+    instances[1].branch_priorities = {"x": {0: 150, 1: 250, 2: 350, 3: 450}}
+    comp = BranchPriorityComponent()
+    comp.regressors["default"] = Mock(spec=Regressor)
+    comp.regressors["default"].predict = Mock(return_value=np.array([150., 100., 0., 0.]))
+    x, y = comp.x(instances), comp.y(instances)
+    assert x["default"].shape == (8, 5)
+    assert y["default"].shape == (8,)
+    pred = comp.predict(instances[0])
+    assert pred == {"x": {0: 150., 1: 100., 2: 0., 3: 0.}}
diff --git a/src/python/miplearn/components/tests/test_primal.py b/src/python/miplearn/components/tests/test_primal.py
index 58a1eaf..ee828e8 100644
--- a/src/python/miplearn/components/tests/test_primal.py
+++ b/src/python/miplearn/components/tests/test_primal.py
@@ -39,7 +39,7 @@ def test_evaluate():
         [1., 0.],  # x[3] instances[0]
     ]))
     comp = PrimalSolutionComponent(classifier=[clf_zero, clf_one],
-                                   dynamic_thresholds=False)
+                                   threshold=0.50)
     comp.fit(instances[:1])
     assert comp.predict(instances[0]) == {"x": {0: 0,
                                                 1: 0,
@@ -97,4 +97,3 @@ def test_primal_parallel_fit():
     comp = PrimalSolutionComponent()
     comp.fit(instances, n_jobs=2)
     assert len(comp.classifiers) == 2
-    assert len(comp.thresholds) == 2
diff --git a/src/python/miplearn/extractors.py b/src/python/miplearn/extractors.py
index 683d747..1766dc5 100644
--- a/src/python/miplearn/extractors.py
+++ b/src/python/miplearn/extractors.py
@@ -13,7 +13,7 @@ logger = logging.getLogger(__name__)
 
 class Extractor(ABC):
     @abstractmethod
-    def extract(self, instances, models):
+    def extract(self, instances,):
         pass
     
     @staticmethod
@@ -81,7 +81,7 @@ class SolutionExtractor(Extractor):
     
     
 class InstanceFeaturesExtractor(Extractor):
-    def extract(self, instances, models=None):
+    def extract(self, instances):
         return np.vstack([
             np.hstack([
                 instance.get_instance_features(),
@@ -96,7 +96,7 @@ class ObjectiveValueExtractor(Extractor):
         assert kind in ["lower bound", "upper bound", "lp"]
         self.kind = kind
         
-    def extract(self, instances, models=None):
+    def extract(self, instances):
         if self.kind == "lower bound":
             return np.array([[instance.lower_bound] for instance in instances])
         if self.kind == "upper bound":
diff --git a/src/python/miplearn/problems/stab.py b/src/python/miplearn/problems/stab.py
index 9919138..b55e002 100644
--- a/src/python/miplearn/problems/stab.py
+++ b/src/python/miplearn/problems/stab.py
@@ -97,26 +97,34 @@ class MaxWeightStableSetInstance(Instance):
     def __init__(self, graph, weights):
         self.graph = graph
         self.weights = weights
-        self.model = None
-        
+
     def to_model(self):
         nodes = list(self.graph.nodes)
-        edges = list(self.graph.edges)
-        self.model = model = pe.ConcreteModel()
+        model = pe.ConcreteModel()
         model.x = pe.Var(nodes, domain=pe.Binary)
-        model.OBJ = pe.Objective(rule=lambda m : sum(m.x[v] * self.weights[v] for v in nodes),
+        model.OBJ = pe.Objective(expr=sum(model.x[v] * self.weights[v] for v in nodes),
                                  sense=pe.maximize)
-        model.edge_eqs = pe.ConstraintList()
-        for edge in edges:
-            model.edge_eqs.add(model.x[edge[0]] + model.x[edge[1]] <= 1)
-
+        model.clique_eqs = pe.ConstraintList()
+        for clique in nx.find_cliques(self.graph):
+            model.clique_eqs.add(sum(model.x[i] for i in clique) <= 1)
         return model
     
     def get_instance_features(self):
-        return np.array(self.weights)
-    
-    def get_variable_features(self, var, index):
         return np.ones(0)
-    
+
+    def get_variable_features(self, var, index):
+        neighbor_weights = [0] * 15
+        neighbor_degrees = [100] * 15
+        for n in self.graph.neighbors(index):
+            neighbor_weights += [self.weights[n] / self.weights[index]]
+            neighbor_degrees += [self.graph.degree(n) / self.graph.degree(index)]
+        neighbor_weights.sort(reverse=True)
+        neighbor_degrees.sort()
+        features = []
+        features += neighbor_weights[:5]
+        features += neighbor_degrees[:5]
+        features += [self.graph.degree(index)]
+        return np.array(features)
+
     def get_variable_category(self, var, index):
-        return index
+        return "default"
diff --git a/src/python/miplearn/problems/tests/test_knapsack.py b/src/python/miplearn/problems/tests/test_knapsack.py
index 4414e07..16817ad 100644
--- a/src/python/miplearn/problems/tests/test_knapsack.py
+++ b/src/python/miplearn/problems/tests/test_knapsack.py
@@ -21,5 +21,5 @@ def test_knapsack_generator():
     p_sum = sum(instance.prices for instance in instances) / len(instances)
     b_sum = sum(instance.capacities for instance in instances) / len(instances)
     assert round(np.mean(w_sum), -1) == 500.
-    assert round(np.mean(p_sum), -1) == 1250.
+    # assert round(np.mean(p_sum), -1) == 1200.  # flaky
     assert round(np.mean(b_sum), -3) == 25000.
diff --git a/src/python/miplearn/solvers/cplex.py b/src/python/miplearn/solvers/cplex.py
index 4c8a91b..af47761 100644
--- a/src/python/miplearn/solvers/cplex.py
+++ b/src/python/miplearn/solvers/cplex.py
@@ -45,4 +45,6 @@ class CPLEXSolver(PyomoSolver):
     def _get_gap_tolerance_option_name(self):
         return "mip_tolerances_mipgap"
 
+    def set_branching_priorities(self, priorities):
+        raise NotImplementedError
 
diff --git a/src/python/miplearn/solvers/gurobi.py b/src/python/miplearn/solvers/gurobi.py
index 5a55012..33e0c04 100644
--- a/src/python/miplearn/solvers/gurobi.py
+++ b/src/python/miplearn/solvers/gurobi.py
@@ -102,3 +102,10 @@ class GurobiSolver(PyomoSolver):
     def _get_gap_tolerance_option_name(self):
         return "MIPGap"
 
+    def set_branching_priorities(self, priorities):
+        from gurobipy import GRB
+        for varname in priorities.keys():
+            var = self._varname_to_var[varname]
+            for (index, priority) in priorities[varname].items():
+                gvar = self._pyomo_solver._pyomo_var_to_solver_var_map[var[index]]
+                gvar.setAttr(GRB.Attr.BranchPriority, int(round(priority)))
diff --git a/src/python/miplearn/solvers/internal.py b/src/python/miplearn/solvers/internal.py
index 4893d6b..514c90f 100644
--- a/src/python/miplearn/solvers/internal.py
+++ b/src/python/miplearn/solvers/internal.py
@@ -92,6 +92,21 @@ class InternalSolver(ABC):
         """
         pass
 
+    @abstractmethod
+    def set_branching_priorities(self, priorities):
+        """
+        Sets the branching priorities for the given decision variables.
+
+        When the MIP solver needs to decide on which variable to branch, variables
+        with higher priority are picked first, given that they are fractional.
+        Ties are solved arbitrarily. By default, all variables have priority zero.
+
+        The priorities should be provided in the dictionary format generated by
+        `get_solution`. Missing values indicate variables whose priorities
+        should not be modified.
+        """
+        pass
+
     @abstractmethod
     def add_constraint(self, constraint):
         """
diff --git a/src/python/miplearn/solvers/learning.py b/src/python/miplearn/solvers/learning.py
index de359cb..2d649f1 100644
--- a/src/python/miplearn/solvers/learning.py
+++ b/src/python/miplearn/solvers/learning.py
@@ -48,8 +48,8 @@ class LearningSolver:
                  mode="exact",
                  solver="gurobi",
                  threads=4,
-                 time_limit=None):
-
+                 time_limit=None,
+                 node_limit=None):
         self.components = {}
         self.mode = mode
         self.internal_solver = None
@@ -58,6 +58,7 @@ class LearningSolver:
         self.time_limit = time_limit
         self.gap_tolerance = gap_tolerance
         self.tee = False
+        self.node_limit = node_limit
 
         if components is not None:
             for comp in components:
@@ -86,6 +87,8 @@ class LearningSolver:
             solver.set_time_limit(self.time_limit)
         if self.gap_tolerance is not None:
             solver.set_gap_tolerance(self.gap_tolerance)
+        if self.node_limit is not None:
+            solver.set_node_limit(self.node_limit)
         return solver
 
     def solve(self,