From e9e34b48d6595ddb70d56e78a986a6212ae6f090 Mon Sep 17 00:00:00 2001 From: Emma Johnson <12833636+emma58@users.noreply.github.com> Date: Sun, 8 Dec 2024 11:55:25 -0500 Subject: [PATCH] Fix hybrid bigm formulation for linear trees (#164) The changes in #163 included changes to the hybrid bigm formulation for linear tree that, while mathematically equivalent, made for a larger formulation in terms of number of constraints. This PR corrects that: It still uses the `gdp.bound_pretransformation` to generate the constraints bounding the features values for each leaf, but it manually transforms the constraints setting the output value to the leaf's linear function, equivalently to @bammari's original implementation. In addition it adds a test to check that the size of the resulting formulation is what is expected. **Legal Acknowledgement**\ By contributing to this software project, I agree my contributions are submitted under the BSD license. I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer. --------- Co-authored-by: Emma Johnson Co-authored-by: jalving --- .../graph_neural_network_formulation.ipynb | 347 +++++++++--------- pyproject.toml | 4 +- src/omlt/__init__.py | 2 +- src/omlt/io/__init__.py | 6 +- src/omlt/linear_tree/lt_formulation.py | 64 +++- src/omlt/neuralnet/__init__.py | 2 +- src/omlt/neuralnet/activations/__init__.py | 8 +- tests/linear_tree/test_lt_formulation.py | 18 +- 8 files changed, 250 insertions(+), 201 deletions(-) diff --git a/docs/notebooks/neuralnet/graph_neural_network_formulation.ipynb b/docs/notebooks/neuralnet/graph_neural_network_formulation.ipynb index 2699d171..e30c94f7 100644 --- a/docs/notebooks/neuralnet/graph_neural_network_formulation.ipynb +++ b/docs/notebooks/neuralnet/graph_neural_network_formulation.ipynb @@ -39,19 +39,9 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 36, "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "2024-05-16 17:32:39.757240: I tensorflow/core/util/port.cc:113] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.\n", - "2024-05-16 17:32:39.808990: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.\n", - "To enable the following instructions: AVX2 AVX512F AVX512_VNNI FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.\n" - ] - } - ], + "outputs": [], "source": [ "import numpy as np\n", "import pyomo.environ as pyo\n", @@ -90,7 +80,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 37, "metadata": {}, "outputs": [ { @@ -135,7 +125,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 38, "metadata": {}, "outputs": [], "source": [ @@ -163,7 +153,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 39, "metadata": {}, "outputs": [ { @@ -171,20 +161,20 @@ "output_type": "stream", "text": [ "Welcome to the CBC MILP Solver \n", - "Version: 2.10.10 \n", - "Build Date: Apr 19 2023 \n", + "Version: 2.10.12 \n", + "Build Date: Sep 3 2024 \n", "\n", - "command line - /opt/conda/bin/cbc -printingOptions all -import /tmp/tmpwsv2x1xb.pyomo.lp -stat=1 -solve -solu /tmp/tmpwsv2x1xb.pyomo.soln (default strategy 1)\n", + "command line - /home/jalving/miniconda3/bin/cbc -printingOptions all -import /tmp/tmplnicse8d.pyomo.lp -stat=1 -solve -solu /tmp/tmplnicse8d.pyomo.soln (default strategy 1)\n", "Option for printingOptions changed from normal to all\n", - "Presolve 172 (-222) rows, 111 (-75) columns and 608 (-267) elements\n", + "Presolve 171 (-223) rows, 109 (-77) columns and 655 (-220) elements\n", "Statistics for presolved model\n", "Original problem has 26 integers (26 of which binary)\n", "Presolved problem has 25 integers (25 of which binary)\n", - "==== 110 zero objective 2 different\n", + "==== 108 zero objective 2 different\n", "1 variables have objective of -0.0421598\n", - "110 variables have objective of 0\n", + "108 variables have objective of 0\n", "==== absolute objective values 2 different\n", - "110 variables have objective of 0\n", + "108 variables have objective of 0\n", "1 variables have objective of 0.0421598\n", "==== for integers 25 zero objective 1 different\n", "25 variables have objective of 0\n", @@ -193,95 +183,95 @@ "===== end objective counts\n", "\n", "\n", - "Problem has 172 rows, 111 columns (1 with objective) and 608 elements\n", + "Problem has 171 rows, 109 columns (1 with objective) and 655 elements\n", "Column breakdown:\n", - "0 of type 0.0->inf, 49 of type 0.0->up, 0 of type lo->inf, \n", - "37 of type lo->up, 0 of type free, 0 of type fixed, \n", + "0 of type 0.0->inf, 44 of type 0.0->up, 0 of type lo->inf, \n", + "40 of type lo->up, 0 of type free, 0 of type fixed, \n", "0 of type -inf->0.0, 0 of type -inf->up, 25 of type 0.0->1.0 \n", "Row breakdown:\n", - "8 of type E 0.0, 0 of type E 1.0, 0 of type E -1.0, \n", - "5 of type E other, 0 of type G 0.0, 0 of type G 1.0, \n", - "0 of type G other, 134 of type L 0.0, 0 of type L 1.0, \n", - "25 of type L other, 0 of type Range 0.0->1.0, 0 of type Range other, \n", + "3 of type E 0.0, 0 of type E 1.0, 0 of type E -1.0, \n", + "8 of type E other, 0 of type G 0.0, 0 of type G 1.0, \n", + "0 of type G other, 133 of type L 0.0, 0 of type L 1.0, \n", + "26 of type L other, 0 of type Range 0.0->1.0, 1 of type Range other, \n", "0 of type Free \n", - "Continuous objective value is 0.315152 - 0.00 seconds\n", + "Continuous objective value is 0.315152 - 0.01 seconds\n", "Cgl0003I 0 fixed, 0 tightened bounds, 2 strengthened rows, 0 substitutions\n", - "Cgl0004I processed model has 166 rows, 105 columns (25 integer (25 of which binary)) and 670 elements\n", - "Cbc0038I Initial state - 5 integers unsatisfied sum - 0.124759\n", - "Cbc0038I Pass 1: suminf. 0.00000 (0) obj. 0.317969 iterations 41\n", + "Cgl0004I processed model has 169 rows, 107 columns (25 integer (25 of which binary)) and 692 elements\n", + "Cbc0038I Initial state - 7 integers unsatisfied sum - 0.178371\n", + "Cbc0038I Pass 1: suminf. 0.00000 (0) obj. 0.317969 iterations 73\n", "Cbc0038I Solution found of 0.317969\n", "Cbc0038I Relaxing continuous gives 0.317969\n", - "Cbc0038I Before mini branch and bound, 20 integers at bound fixed and 63 continuous\n", - "Cbc0038I Full problem 166 rows 105 columns, reduced to 17 rows 13 columns\n", - "Cbc0038I Mini branch and bound did not improve solution (0.02 seconds)\n", - "Cbc0038I Round again with cutoff of 0.317791\n", - "Cbc0038I Pass 2: suminf. 0.00876 (1) obj. 0.317791 iterations 11\n", - "Cbc0038I Pass 3: suminf. 0.18897 (1) obj. 0.317791 iterations 20\n", - "Cbc0038I Pass 4: suminf. 0.00876 (1) obj. 0.317791 iterations 58\n", - "Cbc0038I Pass 5: suminf. 0.18897 (1) obj. 0.317791 iterations 13\n", - "Cbc0038I Pass 6: suminf. 0.00876 (1) obj. 0.317791 iterations 21\n", - "Cbc0038I Pass 7: suminf. 0.00876 (1) obj. 0.317791 iterations 32\n", - "Cbc0038I Pass 8: suminf. 0.18897 (1) obj. 0.317791 iterations 16\n", - "Cbc0038I Pass 9: suminf. 0.00876 (1) obj. 0.317791 iterations 19\n", - "Cbc0038I Pass 10: suminf. 0.00876 (1) obj. 0.317791 iterations 57\n", - "Cbc0038I Pass 11: suminf. 0.18897 (1) obj. 0.317791 iterations 7\n", - "Cbc0038I Pass 12: suminf. 0.00876 (1) obj. 0.317791 iterations 7\n", - "Cbc0038I Pass 13: suminf. 0.00876 (1) obj. 0.317791 iterations 5\n", - "Cbc0038I Pass 14: suminf. 0.18897 (1) obj. 0.317791 iterations 7\n", - "Cbc0038I Pass 15: suminf. 0.00876 (1) obj. 0.317791 iterations 7\n", - "Cbc0038I Pass 16: suminf. 0.00876 (1) obj. 0.317791 iterations 10\n", - "Cbc0038I Pass 17: suminf. 0.18897 (1) obj. 0.317791 iterations 9\n", - "Cbc0038I Pass 18: suminf. 0.00876 (1) obj. 0.317791 iterations 8\n", - "Cbc0038I Pass 19: suminf. 0.00876 (1) obj. 0.317791 iterations 22\n", - "Cbc0038I Pass 20: suminf. 0.18897 (1) obj. 0.317791 iterations 6\n", - "Cbc0038I Pass 21: suminf. 0.00876 (1) obj. 0.317791 iterations 9\n", - "Cbc0038I Pass 22: suminf. 0.00876 (1) obj. 0.317791 iterations 17\n", - "Cbc0038I Pass 23: suminf. 0.18897 (1) obj. 0.317791 iterations 6\n", - "Cbc0038I Pass 24: suminf. 0.00876 (1) obj. 0.317791 iterations 5\n", - "Cbc0038I Pass 25: suminf. 0.00876 (1) obj. 0.317791 iterations 10\n", - "Cbc0038I Pass 26: suminf. 0.18897 (1) obj. 0.317791 iterations 6\n", - "Cbc0038I Pass 27: suminf. 0.00876 (1) obj. 0.317791 iterations 5\n", - "Cbc0038I Pass 28: suminf. 0.00876 (1) obj. 0.317791 iterations 30\n", - "Cbc0038I Pass 29: suminf. 0.18897 (1) obj. 0.317791 iterations 5\n", - "Cbc0038I Pass 30: suminf. 0.00876 (1) obj. 0.317791 iterations 6\n", - "Cbc0038I Pass 31: suminf. 0.00876 (1) obj. 0.317791 iterations 3\n", + "Cbc0038I Before mini branch and bound, 18 integers at bound fixed and 59 continuous\n", + "Cbc0038I Full problem 169 rows 107 columns, reduced to 33 rows 17 columns\n", + "Cbc0038I Mini branch and bound did not improve solution (0.33 seconds)\n", + "Cbc0038I Round again with cutoff of 0.317678\n", + "Cbc0038I Pass 2: suminf. 0.01308 (1) obj. 0.317678 iterations 20\n", + "Cbc0038I Pass 3: suminf. 0.28789 (1) obj. 0.317678 iterations 7\n", + "Cbc0038I Pass 4: suminf. 0.01308 (1) obj. 0.317678 iterations 23\n", + "Cbc0038I Pass 5: suminf. 0.28789 (1) obj. 0.317678 iterations 10\n", + "Cbc0038I Pass 6: suminf. 0.01308 (1) obj. 0.317678 iterations 8\n", + "Cbc0038I Pass 7: suminf. 0.01308 (1) obj. 0.317678 iterations 30\n", + "Cbc0038I Pass 8: suminf. 0.28789 (1) obj. 0.317678 iterations 8\n", + "Cbc0038I Pass 9: suminf. 0.01308 (1) obj. 0.317678 iterations 6\n", + "Cbc0038I Pass 10: suminf. 0.01308 (1) obj. 0.317678 iterations 16\n", + "Cbc0038I Pass 11: suminf. 0.28789 (1) obj. 0.317678 iterations 3\n", + "Cbc0038I Pass 12: suminf. 0.01308 (1) obj. 0.317678 iterations 2\n", + "Cbc0038I Pass 13: suminf. 0.01308 (1) obj. 0.317678 iterations 21\n", + "Cbc0038I Pass 14: suminf. 0.28789 (1) obj. 0.317678 iterations 6\n", + "Cbc0038I Pass 15: suminf. 0.01308 (1) obj. 0.317678 iterations 4\n", + "Cbc0038I Pass 16: suminf. 0.01308 (1) obj. 0.317678 iterations 9\n", + "Cbc0038I Pass 17: suminf. 0.28789 (1) obj. 0.317678 iterations 6\n", + "Cbc0038I Pass 18: suminf. 0.01308 (1) obj. 0.317678 iterations 9\n", + "Cbc0038I Pass 19: suminf. 0.01308 (1) obj. 0.317678 iterations 7\n", + "Cbc0038I Pass 20: suminf. 0.28789 (1) obj. 0.317678 iterations 3\n", + "Cbc0038I Pass 21: suminf. 0.01308 (1) obj. 0.317678 iterations 2\n", + "Cbc0038I Pass 22: suminf. 0.01308 (1) obj. 0.317678 iterations 16\n", + "Cbc0038I Pass 23: suminf. 0.28789 (1) obj. 0.317678 iterations 2\n", + "Cbc0038I Pass 24: suminf. 0.01308 (1) obj. 0.317678 iterations 1\n", + "Cbc0038I Pass 25: suminf. 0.01308 (1) obj. 0.317678 iterations 29\n", + "Cbc0038I Pass 26: suminf. 0.28789 (1) obj. 0.317678 iterations 2\n", + "Cbc0038I Pass 27: suminf. 0.01308 (1) obj. 0.317678 iterations 1\n", + "Cbc0038I Pass 28: suminf. 0.01308 (1) obj. 0.317678 iterations 10\n", + "Cbc0038I Pass 29: suminf. 0.28789 (1) obj. 0.317678 iterations 4\n", + "Cbc0038I Pass 30: suminf. 0.01308 (1) obj. 0.317678 iterations 3\n", + "Cbc0038I Pass 31: suminf. 0.01308 (1) obj. 0.317678 iterations 3\n", "Cbc0038I No solution found this major pass\n", - "Cbc0038I Before mini branch and bound, 1 integers at bound fixed and 47 continuous\n", - "Cbc0038I Full problem 166 rows 105 columns, reduced to 48 rows 27 columns\n", - "Cbc0038I Mini branch and bound did not improve solution (0.04 seconds)\n", - "Cbc0038I After 0.04 seconds - Feasibility pump exiting with objective of 0.317969 - took 0.03 seconds\n", - "Cbc0012I Integer solution of 0.31796885 found by feasibility pump after 0 iterations and 0 nodes (0.05 seconds)\n", - "Cbc0038I Full problem 166 rows 105 columns, reduced to 48 rows 27 columns\n", - "Cbc0031I 3 added rows had average density of 3.3333333\n", - "Cbc0013I At root node, 31 cuts changed objective from 0.31628066 to 0.31796885 in 1 passes\n", - "Cbc0014I Cut generator 0 (Probing) - 19 row cuts average 3.0 elements, 1 column cuts (1 active) in 0.000 seconds - new frequency is 1\n", - "Cbc0014I Cut generator 1 (Gomory) - 3 row cuts average 8.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is 1\n", + "Cbc0038I Before mini branch and bound, 2 integers at bound fixed and 49 continuous\n", + "Cbc0038I Full problem 169 rows 107 columns, reduced to 56 rows 27 columns\n", + "Cbc0038I Mini branch and bound did not improve solution (0.65 seconds)\n", + "Cbc0038I After 0.65 seconds - Feasibility pump exiting with objective of 0.317969 - took 0.44 seconds\n", + "Cbc0012I Integer solution of 0.31796885 found by feasibility pump after 0 iterations and 0 nodes (0.66 seconds)\n", + "Cbc0038I Full problem 169 rows 107 columns, reduced to 56 rows 27 columns\n", + "Cbc0031I 4 added rows had average density of 8.5\n", + "Cbc0013I At root node, 24 cuts changed objective from 0.31515217 to 0.31796885 in 1 passes\n", + "Cbc0014I Cut generator 0 (Probing) - 14 row cuts average 3.0 elements, 1 column cuts (1 active) in 0.000 seconds - new frequency is 1\n", + "Cbc0014I Cut generator 1 (Gomory) - 3 row cuts average 10.3 elements, 0 column cuts (0 active) in 0.001 seconds - new frequency is 1\n", "Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100\n", "Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100\n", - "Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 3 row cuts average 3.3 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is 1\n", + "Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 2 row cuts average 3.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is 1\n", "Cbc0014I Cut generator 5 (FlowCover) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100\n", - "Cbc0014I Cut generator 6 (TwoMirCuts) - 6 row cuts average 6.2 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is 1\n", - "Cbc0001I Search completed - best objective 0.3179688539269278, took 17 iterations and 0 nodes (0.06 seconds)\n", + "Cbc0014I Cut generator 6 (TwoMirCuts) - 5 row cuts average 7.8 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is 1\n", + "Cbc0001I Search completed - best objective 0.3179688462393472, took 20 iterations and 0 nodes (0.75 seconds)\n", "Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost\n", - "Cuts at root node changed objective from 0.316281 to 0.317969\n", - "Probing was tried 1 times and created 20 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", - "Gomory was tried 1 times and created 3 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", + "Cuts at root node changed objective from 0.315152 to 0.317969\n", + "Probing was tried 1 times and created 15 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", + "Gomory was tried 1 times and created 3 cuts of which 0 were active after adding rounds of cuts (0.001 seconds)\n", "Knapsack was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", "Clique was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", - "MixedIntegerRounding2 was tried 1 times and created 3 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", + "MixedIntegerRounding2 was tried 1 times and created 2 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", "FlowCover was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", - "TwoMirCuts was tried 1 times and created 6 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", + "TwoMirCuts was tried 1 times and created 5 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", "ZeroHalf was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", "\n", "Result - Optimal solution found\n", "\n", "Objective value: 0.31796885\n", "Enumerated nodes: 0\n", - "Total iterations: 17\n", - "Time (CPU seconds): 0.06\n", - "Time (Wallclock seconds): 0.03\n", + "Total iterations: 20\n", + "Time (CPU seconds): 0.84\n", + "Time (Wallclock seconds): 0.07\n", "\n", - "Total time (CPU seconds): 0.07 (Wallclock seconds): 0.04\n", + "Total time (CPU seconds): 0.93 (Wallclock seconds): 0.08\n", "\n" ] } @@ -312,15 +302,28 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 32, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "[[0.31796885]]\n" - ] + "data": { + "text/plain": [ + "Sequential(\n", + " (0) - GCNConv(2, 4): x, edge_index -> x\n", + " (1) - ReLU(): x -> x\n", + " (2) - GCNConv(4, 4): x, edge_index -> x\n", + " (3) - ReLU(): x -> x\n", + " (4) - Linear(in_features=4, out_features=4, bias=True): x -> x\n", + " (5) - : x, None -> x\n", + " (6) - Linear(in_features=4, out_features=2, bias=True): x -> x\n", + " (7) - ReLU(): x -> x\n", + " (8) - Linear(in_features=2, out_features=1, bias=True): x -> x\n", + ")" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ @@ -335,7 +338,9 @@ "X = np.array(X).reshape(3, 2)\n", "edges = np.transpose(np.array(edges)).reshape(2, -1)\n", "nn.eval()\n", - "print(nn1(torch.tensor(X).float(), torch.tensor(edges)).detach().numpy())" + "# TODO @zshiqiang: update for latest torch.geometric # noqa: FIX002\n", + "# https://github.com/cog-imperial/OMLT/issues/166\n", + "# print(nn1(torch.tensor(X).float(), torch.tensor(edges)).detach().numpy()) : noqa: ERA001" ] }, { @@ -347,7 +352,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 40, "metadata": {}, "outputs": [ { @@ -362,7 +367,7 @@ " For more information visit https://github.com/coin-or/Ipopt\n", "******************************************************************************\n", "\n", - "This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.1.\n", + "This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.\n", "\n", "Number of nonzeros in equality constraint Jacobian...: 395\n", "Number of nonzeros in inequality constraint Jacobian.: 276\n", @@ -380,65 +385,61 @@ "\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 0 0.0000000e+00 4.73e-01 1.32e-04 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", - " 1 7.5562415e-02 3.99e-01 6.07e+01 -1.0 6.41e-01 - 4.88e-02 1.57e-01f 1\n", - " 2 1.6790548e-01 3.08e-01 4.25e+01 -1.0 4.05e-01 - 1.81e-01 2.28e-01f 1\n", - " 3 2.5794319e-01 2.19e-01 3.18e+01 -1.0 3.13e-01 - 5.10e-01 2.88e-01f 1\n", - " 4 3.9100053e-01 8.85e-02 1.45e+01 -1.0 2.23e-01 - 9.76e-01 5.97e-01h 1\n", - " 5 4.4307883e-01 3.72e-02 2.75e+01 -1.0 1.41e-01 - 1.00e+00 5.79e-01h 1\n", - " 6 4.6540208e-01 1.52e-02 6.33e+01 -1.0 8.16e-02 - 1.00e+00 5.91e-01h 1\n", + " 1 7.5562414e-02 3.99e-01 6.07e+01 -1.0 6.41e-01 - 4.88e-02 1.57e-01f 1\n", + " 2 1.6790546e-01 3.08e-01 4.25e+01 -1.0 4.05e-01 - 1.81e-01 2.28e-01f 1\n", + " 3 2.5794318e-01 2.19e-01 3.18e+01 -1.0 3.13e-01 - 5.10e-01 2.88e-01f 1\n", + " 4 3.9100052e-01 8.85e-02 1.45e+01 -1.0 2.23e-01 - 9.76e-01 5.97e-01h 1\n", + " 5 4.4307882e-01 3.72e-02 2.75e+01 -1.0 1.41e-01 - 1.00e+00 5.79e-01h 1\n", + " 6 4.6540207e-01 1.52e-02 6.33e+01 -1.0 8.16e-02 - 1.00e+00 5.91e-01h 1\n", " 7 4.7445376e-01 6.31e-03 1.55e+02 -1.0 3.67e-02 - 1.00e+00 5.85e-01h 1\n", " 8 4.7820844e-01 2.61e-03 3.74e+02 -1.0 1.47e-02 - 1.00e+00 5.86e-01h 1\n", - " 9 4.7976341e-01 1.08e-03 9.04e+02 -1.0 6.36e-03 - 1.00e+00 5.86e-01h 1\n", + " 9 4.7976340e-01 1.08e-03 9.04e+02 -1.0 6.36e-03 - 1.00e+00 5.86e-01h 1\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", " 10 4.8040743e-01 4.48e-04 2.18e+03 -1.0 2.53e-03 - 1.00e+00 5.86e-01h 1\n", - " 11 4.8067424e-01 1.86e-04 5.26e+03 -1.0 1.09e-03 - 1.00e+00 5.86e-01h 1\n", + " 11 4.8067423e-01 1.86e-04 5.26e+03 -1.0 1.09e-03 - 1.00e+00 5.86e-01h 1\n", " 12 4.8078473e-01 7.67e-05 1.27e+04 -1.0 4.34e-04 - 1.00e+00 5.87e-01h 1\n", " 13 4.8083051e-01 3.16e-05 3.04e+04 -1.0 1.87e-04 - 1.00e+00 5.88e-01h 1\n", - " 14 4.8084947e-01 1.29e-05 7.21e+04 -1.0 7.40e-05 - 1.00e+00 5.91e-01h 1\n", + " 14 4.8084946e-01 1.29e-05 7.21e+04 -1.0 7.40e-05 - 1.00e+00 5.91e-01h 1\n", " 15 4.8085732e-01 5.18e-06 1.67e+05 -1.0 3.15e-05 - 1.00e+00 5.99e-01h 1\n", " 16 4.8086057e-01 1.98e-06 3.65e+05 -1.0 1.21e-05 - 1.00e+00 6.18e-01h 1\n", - " 17 4.8086191e-01 6.64e-07 6.79e+05 -1.0 4.84e-06 - 1.00e+00 6.65e-01h 1\n", + " 17 4.8086190e-01 6.64e-07 6.79e+05 -1.0 4.84e-06 - 1.00e+00 6.65e-01h 1\n", " 18 4.8086192e-01 6.47e-07 3.49e+06 -1.0 1.54e-06 - 1.00e+00 2.47e-02f 6\n", - " 19 4.8086258e-01 3.98e-10 1.00e-06 -1.0 1.52e-06 - 1.00e+00 1.00e+00h 1\n", + " 19 4.8086258e-01 1.92e-10 1.00e-06 -1.0 1.52e-06 - 1.00e+00 1.00e+00h 1\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", - " 20 4.8086253e-01 2.78e-09 4.52e+02 -8.6 7.22e-05 - 1.00e+00 1.00e+00h 1\n", - " 21 4.8001912e-01 2.86e-02 2.80e+02 -8.6 1.17e+00 - 5.23e-01 1.00e+00f 1\n", + " 20 4.8086252e-01 1.96e-09 4.52e+02 -8.6 7.22e-05 - 1.00e+00 1.00e+00h 1\n", + " 21 4.8001911e-01 2.86e-02 2.80e+02 -8.6 1.17e+00 - 5.23e-01 1.00e+00f 1\n", " 22 4.8001743e-01 1.08e-02 2.14e+01 -8.6 2.15e-01 - 9.00e-01 6.57e-01h 1\n", - " 23 4.8001271e-01 1.79e-03 2.27e+01 -8.6 3.03e-01 - 8.31e-01 1.00e+00h 1\n", - " 24 4.8000768e-01 1.81e-04 4.21e+01 -8.6 9.74e-02 - 7.43e-01 1.00e+00h 1\n", - " 25 4.8000767e-01 1.79e-04 5.50e+02 -8.6 1.68e-02 - 1.00e+00 9.31e-03h 1\n", - " 26 4.8000669e-01 5.00e-06 9.43e-09 -8.6 1.76e-02 - 1.00e+00 1.00e+00f 1\n", - " 27 4.8000670e-01 7.87e-11 8.32e+01 -8.6 7.60e-05 - 4.17e-01 1.00e+00h 1\n", - " 28 4.8000670e-01 7.86e-11 1.41e+02 -8.6 8.40e-06 - 1.00e+00 1.93e-04h 1\n", - " 29 4.8000670e-01 3.40e-10 1.73e+02 -8.6 7.68e-06 - 6.59e-02 1.00e+00f 1\n", + " 23 4.8001271e-01 1.79e-03 2.28e+01 -8.6 3.03e-01 - 8.31e-01 1.00e+00h 1\n", + " 24 4.8000768e-01 1.81e-04 4.22e+01 -8.6 9.74e-02 - 7.43e-01 1.00e+00h 1\n", + " 25 4.8000768e-01 1.81e-04 3.59e+01 -8.6 9.42e-03 -4.0 4.08e-01 4.08e-04h 1\n", + " 26 4.8000767e-01 1.80e-04 8.64e+01 -8.6 2.02e-02 - 4.83e-01 2.60e-03f 2\n", + " 27 4.8000767e-01 1.80e-04 1.41e+02 -8.6 1.89e-02 - 1.00e+00 5.52e-04h 1\n", + " 28 4.8000668e-01 5.07e-06 5.74e+00 -8.6 1.78e-02 - 2.96e-01 1.00e+00f 1\n", + " 29 4.8000668e-01 5.06e-06 1.41e+02 -8.6 3.89e-05 - 1.00e+00 9.95e-04h 1\n", "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", - " 30 4.8000670e-01 3.40e-10 1.74e+02 -8.6 2.35e-05 - 1.00e+00 2.85e-04h 1\n", - " 31 4.8000670e-01 5.51e-09 1.32e+00 -8.6 3.08e-08 - 2.22e-01 1.00e+00f 1\n", - " 32 4.8000670e-01 9.70e-11 1.23e+02 -8.6 1.11e-04 - 8.72e-01 1.82e-04h 2\n", - " 33 4.8000670e-01 1.52e-10 1.76e-10 -8.6 1.50e-08 - 1.00e+00 1.00e+00h 1\n", + " 30 4.8000669e-01 9.53e-11 1.99e-10 -8.6 4.54e-05 - 1.00e+00 1.00e+00f 1\n", "\n", - "Number of Iterations....: 33\n", + "Number of Iterations....: 30\n", "\n", " (scaled) (unscaled)\n", - "Objective...............: 4.8000669509919508e-01 4.8000669509919508e-01\n", - "Dual infeasibility......: 1.7627861836399183e-10 1.7627861836399183e-10\n", - "Constraint violation....: 1.5176976342345938e-10 1.5176976342345938e-10\n", + "Objective...............: 4.8000669067803564e-01 4.8000669067803564e-01\n", + "Dual infeasibility......: 1.9944890902928439e-10 1.9944890902928439e-10\n", + "Constraint violation....: 9.5274010902812734e-11 9.5274010902812734e-11\n", "Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00\n", - "Complementarity.........: 3.0531777979505568e-09 3.0531777979505568e-09\n", - "Overall NLP error.......: 1.5176976342345938e-10 3.0531777979505568e-09\n", + "Complementarity.........: 2.5481058108255804e-09 2.5481058108255804e-09\n", + "Overall NLP error.......: 9.5274010902812734e-11 2.5481058108255804e-09\n", "\n", "\n", - "Number of objective function evaluations = 40\n", - "Number of objective gradient evaluations = 34\n", - "Number of equality constraint evaluations = 40\n", - "Number of inequality constraint evaluations = 40\n", - "Number of equality constraint Jacobian evaluations = 34\n", - "Number of inequality constraint Jacobian evaluations = 34\n", - "Number of Lagrangian Hessian evaluations = 33\n", - "Total seconds in IPOPT = 0.052\n", + "Number of objective function evaluations = 37\n", + "Number of objective gradient evaluations = 31\n", + "Number of equality constraint evaluations = 37\n", + "Number of inequality constraint evaluations = 37\n", + "Number of equality constraint Jacobian evaluations = 31\n", + "Number of inequality constraint Jacobian evaluations = 31\n", + "Number of Lagrangian Hessian evaluations = 30\n", + "Total seconds in IPOPT = 0.059\n", "\n", - "EXIT: Optimal Solution Found.\n", - "\b" + "EXIT: Optimal Solution Found.\n" ] } ], @@ -473,7 +474,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 41, "metadata": {}, "outputs": [], "source": [ @@ -514,7 +515,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 42, "metadata": {}, "outputs": [ { @@ -522,21 +523,21 @@ "output_type": "stream", "text": [ "Welcome to the CBC MILP Solver \n", - "Version: 2.10.10 \n", - "Build Date: Apr 19 2023 \n", + "Version: 2.10.12 \n", + "Build Date: Sep 3 2024 \n", "\n", - "command line - /opt/conda/bin/cbc -printingOptions all -import /tmp/tmp0s5lbbp6.pyomo.lp -stat=1 -solve -solu /tmp/tmp0s5lbbp6.pyomo.soln (default strategy 1)\n", + "command line - /home/jalving/miniconda3/bin/cbc -printingOptions all -import /tmp/tmpe02s_bpw.pyomo.lp -stat=1 -solve -solu /tmp/tmpe02s_bpw.pyomo.soln (default strategy 1)\n", "Option for printingOptions changed from normal to all\n", - "Presolve 260 (-137) rows, 141 (-51) columns and 852 (-173) elements\n", + "Presolve 266 (-131) rows, 147 (-45) columns and 875 (-150) elements\n", "Statistics for presolved model\n", "Original problem has 32 integers (32 of which binary)\n", "Presolved problem has 29 integers (29 of which binary)\n", - "==== 139 zero objective 3 different\n", - "139 variables have objective of 0\n", + "==== 145 zero objective 3 different\n", + "145 variables have objective of 0\n", "1 variables have objective of 0.203177\n", "1 variables have objective of 0.686721\n", "==== absolute objective values 3 different\n", - "139 variables have objective of 0\n", + "145 variables have objective of 0\n", "1 variables have objective of 0.203177\n", "1 variables have objective of 0.686721\n", "==== for integers 29 zero objective 1 different\n", @@ -546,28 +547,29 @@ "===== end objective counts\n", "\n", "\n", - "Problem has 260 rows, 141 columns (2 with objective) and 852 elements\n", + "Problem has 266 rows, 147 columns (2 with objective) and 875 elements\n", "Column breakdown:\n", - "0 of type 0.0->inf, 62 of type 0.0->up, 0 of type lo->inf, \n", - "50 of type lo->up, 0 of type free, 0 of type fixed, \n", + "0 of type 0.0->inf, 61 of type 0.0->up, 0 of type lo->inf, \n", + "57 of type lo->up, 0 of type free, 0 of type fixed, \n", "0 of type -inf->0.0, 0 of type -inf->up, 29 of type 0.0->1.0 \n", "Row breakdown:\n", - "0 of type E 0.0, 0 of type E 1.0, 0 of type E -1.0, \n", - "26 of type E other, 0 of type G 0.0, 0 of type G 1.0, \n", - "0 of type G other, 154 of type L 0.0, 24 of type L 1.0, \n", - "56 of type L other, 0 of type Range 0.0->1.0, 0 of type Range other, \n", + "2 of type E 0.0, 0 of type E 1.0, 0 of type E -1.0, \n", + "30 of type E other, 0 of type G 0.0, 0 of type G 1.0, \n", + "0 of type G other, 152 of type L 0.0, 24 of type L 1.0, \n", + "58 of type L other, 0 of type Range 0.0->1.0, 0 of type Range other, \n", "0 of type Free \n", - "Continuous objective value is 0.107106 - 0.00 seconds\n", - "Cgl0004I processed model has 237 rows, 118 columns (29 integer (29 of which binary)) and 969 elements\n", - "Cbc0038I Initial state - 17 integers unsatisfied sum - 1.61761\n", - "Cbc0038I Pass 1: suminf. 1.01765 (9) obj. 0.107106 iterations 46\n", + "Continuous objective value is 0.107106 - 0.03 seconds\n", + "Cgl0003I 0 fixed, 0 tightened bounds, 4 strengthened rows, 0 substitutions\n", + "Cgl0004I processed model has 250 rows, 131 columns (29 integer (29 of which binary)) and 985 elements\n", + "Cbc0038I Initial state - 19 integers unsatisfied sum - 2.07755\n", + "Cbc0038I Pass 1: suminf. 1.01765 (9) obj. 0.107106 iterations 53\n", "Cbc0038I Solution found of 0.107106\n", "Cbc0038I Relaxing continuous gives 0.107106\n", - "Cbc0038I Before mini branch and bound, 12 integers at bound fixed and 38 continuous\n", - "Cbc0038I Mini branch and bound did not improve solution (0.01 seconds)\n", - "Cbc0038I After 0.02 seconds - Feasibility pump exiting with objective of 0.107106 - took 0.01 seconds\n", - "Cbc0012I Integer solution of 0.10710584 found by feasibility pump after 0 iterations and 0 nodes (0.02 seconds)\n", - "Cbc0001I Search completed - best objective 0.1071058437228203, took 0 iterations and 0 nodes (0.02 seconds)\n", + "Cbc0038I Before mini branch and bound, 9 integers at bound fixed and 37 continuous\n", + "Cbc0038I Mini branch and bound did not improve solution (0.11 seconds)\n", + "Cbc0038I After 0.13 seconds - Feasibility pump exiting with objective of 0.107106 - took 0.05 seconds\n", + "Cbc0012I Integer solution of 0.10710584 found by feasibility pump after 0 iterations and 0 nodes (0.13 seconds)\n", + "Cbc0001I Search completed - best objective 0.10710584, took 0 iterations and 0 nodes (0.15 seconds)\n", "Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost\n", "Cuts at root node changed objective from 0.107106 to 0.107106\n", "Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", @@ -584,10 +586,10 @@ "Objective value: 0.10710584\n", "Enumerated nodes: 0\n", "Total iterations: 0\n", - "Time (CPU seconds): 0.03\n", - "Time (Wallclock seconds): 0.02\n", + "Time (CPU seconds): 0.20\n", + "Time (Wallclock seconds): 0.01\n", "\n", - "Total time (CPU seconds): 0.03 (Wallclock seconds): 0.02\n", + "Total time (CPU seconds): 0.23 (Wallclock seconds): 0.02\n", "\n" ] } @@ -620,18 +622,11 @@ "# solve the optimization problem\n", "status = pyo.SolverFactory(\"cbc\").solve(m3, tee=True)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -645,9 +640,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.13" + "version": "3.12.3" } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } diff --git a/pyproject.toml b/pyproject.toml index 2b215f71..b10b1a84 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -74,12 +74,12 @@ extend-exclude = ["src/omlt/_version.py"] [tool.ruff.lint] select = ["ALL"] ignore = [ - "ANN101", "ANN401", "COM812", "ISC001", "SLF001", "ARG001", + "PLC0206", "N803", "N806", # Remove these after issue https://github.com/cog-imperial/OMLT/issues/153 is fixed. @@ -96,7 +96,7 @@ ignore = [ "ANN002", "ANN201", "ANN202", - "ANN204", + "ANN204" ] [tool.ruff.lint.pydocstyle] diff --git a/src/omlt/__init__.py b/src/omlt/__init__.py index 30fe7b9c..22daa160 100644 --- a/src/omlt/__init__.py +++ b/src/omlt/__init__.py @@ -13,7 +13,7 @@ from omlt.scaling import OffsetScaling __all__ = [ - "OmltBlock", "OffsetScaling", + "OmltBlock", "__version__", ] diff --git a/src/omlt/io/__init__.py b/src/omlt/io/__init__.py index 64fa72e1..9a1fe7ee 100644 --- a/src/omlt/io/__init__.py +++ b/src/omlt/io/__init__.py @@ -17,11 +17,11 @@ __all__ = [ "keras_available", + "load_keras_sequential", + "load_onnx_neural_network", + "load_onnx_neural_network_with_bounds", "onnx_available", "torch_available", "torch_geometric_available", - "load_onnx_neural_network", - "load_onnx_neural_network_with_bounds", "write_onnx_model_with_bounds", - "load_keras_sequential", ] diff --git a/src/omlt/linear_tree/lt_formulation.py b/src/omlt/linear_tree/lt_formulation.py index d72df160..581719e4 100644 --- a/src/omlt/linear_tree/lt_formulation.py +++ b/src/omlt/linear_tree/lt_formulation.py @@ -106,6 +106,7 @@ def _build_formulation(self): output_vars=self.block.scaled_outputs, transformation=self.transformation, epsilon=self.epsilon, + include_leaf_equalities=True, ) @@ -170,23 +171,51 @@ def _build_formulation(self): This method is called by the OmltBlock to build the corresponding mathematical formulation on the Pyomo block. """ + block = self.block + leaves = self.model_definition.leaves + _setup_scaled_inputs_outputs( - self.block, + block, self.model_definition.scaling_object, self.model_definition.scaled_input_bounds, ) + input_vars = self.block.scaled_inputs + _add_gdp_formulation_to_block( - block=self.block, + block=block, model_definition=self.model_definition, - input_vars=self.block.scaled_inputs, + input_vars=input_vars, output_vars=self.block.scaled_outputs, transformation="custom", epsilon=self.epsilon, + include_leaf_equalities=False, ) - pe.TransformationFactory("gdp.bound_pretransformation").apply_to(self.block) - pe.TransformationFactory("gdp.binary_multiplication").apply_to(self.block) + pe.TransformationFactory("gdp.bound_pretransformation").apply_to(block) + # It doesn't really matter what transformation we call next, so we just + # use bigm--all it's going to do is create the exactly-one constraints + # and mark all the disjunctive parts of the model as transformed. + pe.TransformationFactory("gdp.bigm").apply_to(block) + + # We now create the \sum((a_l^Tx + b_l)*y_l for l in leaves) = d constraints + # manually. + features = np.arange(0, self.model_definition.n_inputs) + + @block.Constraint(list(leaves.keys())) + def linear_constraint(mdl, tree): + leaf_ids = list(leaves[tree].keys()) + return block.intermediate_output[tree] == sum( + ( + sum( + leaves[tree][leaf]["slope"][feat] * input_vars[feat] + for feat in features + ) + + leaves[tree][leaf]["intercept"] + ) + * block.disjunct[tree, leaf].binary_indicator_var + for leaf in leaf_ids + ) def _build_output_bounds(model_def, input_bounds): @@ -232,7 +261,13 @@ def _build_output_bounds(model_def, input_bounds): def _add_gdp_formulation_to_block( # noqa: PLR0913 - block, model_definition, input_vars, output_vars, transformation, epsilon + block, + model_definition, + input_vars, + output_vars, + transformation, + epsilon, + include_leaf_equalities, ): """This function adds the GDP representation to the OmltBlock using Pyomo.GDP. @@ -245,7 +280,9 @@ def _add_gdp_formulation_to_block( # noqa: PLR0913 epsilon: Tolerance to use in enforcing that choosing the right branch of a linear tree node can only happen if the feature is strictly greater than the branch value. - + include_leaf_equalities: boolean to indicate if the formulation + should include the equalities setting the leaf values or not. + (default: True) """ leaves = model_definition.leaves input_bounds = model_definition.scaled_input_bounds @@ -283,12 +320,13 @@ def ub_rule(dsj, feat): dsj.ub_constraint = pe.Constraint(features, rule=ub_rule) - slope = leaves[tree][leaf]["slope"] - intercept = leaves[tree][leaf]["intercept"] - dsj.linear_exp = pe.Constraint( - expr=sum(slope[k] * input_vars[k] for k in features) + intercept - == block.intermediate_output[tree] - ) + if include_leaf_equalities: + slope = leaves[tree][leaf]["slope"] + intercept = leaves[tree][leaf]["intercept"] + dsj.linear_exp = pe.Constraint( + expr=sum(slope[k] * input_vars[k] for k in features) + intercept + == block.intermediate_output[tree] + ) block.disjunct = Disjunct(t_l, rule=disjuncts_rule) diff --git a/src/omlt/neuralnet/__init__.py b/src/omlt/neuralnet/__init__.py index 014de739..e4072c30 100644 --- a/src/omlt/neuralnet/__init__.py +++ b/src/omlt/neuralnet/__init__.py @@ -32,9 +32,9 @@ ) __all__ = [ - "NetworkDefinition", "FullSpaceNNFormulation", "FullSpaceSmoothNNFormulation", + "NetworkDefinition", "ReducedSpaceNNFormulation", "ReducedSpaceSmoothNNFormulation", "ReluBigMFormulation", diff --git a/src/omlt/neuralnet/activations/__init__.py b/src/omlt/neuralnet/activations/__init__.py index d60d00c4..c3198e54 100644 --- a/src/omlt/neuralnet/activations/__init__.py +++ b/src/omlt/neuralnet/activations/__init__.py @@ -30,16 +30,16 @@ NON_INCREASING_ACTIVATIONS: list[Any] = [] __all__ = [ - "linear_activation_constraint", - "linear_activation_function", + "ACTIVATION_FUNCTION_MAP", + "NON_INCREASING_ACTIVATIONS", "ComplementarityReLUActivation", "bigm_relu_activation_constraint", + "linear_activation_constraint", + "linear_activation_function", "sigmoid_activation_constraint", "sigmoid_activation_function", "softplus_activation_constraint", "softplus_activation_function", "tanh_activation_constraint", "tanh_activation_function", - "ACTIVATION_FUNCTION_MAP", - "NON_INCREASING_ACTIVATIONS", ] diff --git a/tests/linear_tree/test_lt_formulation.py b/tests/linear_tree/test_lt_formulation.py index db9121dd..1053c8a0 100644 --- a/tests/linear_tree/test_lt_formulation.py +++ b/tests/linear_tree/test_lt_formulation.py @@ -1,6 +1,8 @@ import numpy as np import pyomo.environ as pe import pytest +from pyomo.common.collections import ComponentSet +from pyomo.core.expr import identify_variables from omlt.dependencies import lineartree_available @@ -245,7 +247,7 @@ def test_nonzero_epsilon(): solution = (pe.value(model_good.x), pe.value(model_good.y)) y_pred = regr_small.predict(np.array(solution[0]).reshape(1, -1)) # With epsilon, the model matches the tree prediction - assert y_pred[0] == pytest.approx(solution[1]) + assert y_pred[0] == pytest.approx(solution[1], abs=1e-4) @pytest.mark.skipif( @@ -657,6 +659,20 @@ def test_hybrid_bigm_formulation_multi_var(): model1.lt = OmltBlock() model1.lt.build_formulation(formulation1_lt) + num_constraints = 0 + var_set = ComponentSet() + for cons in model1.lt.component_data_objects(pe.Constraint, active=True): + num_constraints += 1 + for v in identify_variables(cons.expr): + var_set.add(v) + + num_leaves = len(ltmodel_small.leaves[0]) + # binary for each leaf + two inputs and an output + 5 scaled input/output vars + assert len(var_set) == num_leaves + 3 + 4 + # 2 bounds constraints for each input, the xor, the output constraint, and + # four scaling constraints from OMLT + assert num_constraints == 2 * 2 + 1 + 1 + 4 + @model1.Constraint() def connect_input1(mdl): return mdl.x0 == mdl.lt.inputs[0]