1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the program and library             */
4    	/*         SCIP --- Solving Constraint Integer Programs                      */
5    	/*                                                                           */
6    	/*  Copyright 2002-2023 Zuse Institute Berlin                                */
7    	/*                                                                           */
8    	/*  Licensed under the Apache License, Version 2.0 (the "License");          */
9    	/*  you may not use this file except in compliance with the License.         */
10   	/*  You may obtain a copy of the License at                                  */
11   	/*                                                                           */
12   	/*      http://www.apache.org/licenses/LICENSE-2.0                           */
13   	/*                                                                           */
14   	/*  Unless required by applicable law or agreed to in writing, software      */
15   	/*  distributed under the License is distributed on an "AS IS" BASIS,        */
16   	/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17   	/*  See the License for the specific language governing permissions and      */
18   	/*  limitations under the License.                                           */
19   	/*                                                                           */
20   	/*  You should have received a copy of the Apache-2.0 license                */
21   	/*  along with SCIP; see the file LICENSE. If not visit scipopt.org.         */
22   	/*                                                                           */
23   	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24   	
25   	/**@file   branch_gomory.c
26   	 * @ingroup DEFPLUGINS_BRANCH
27   	 * @brief  Gomory cut branching rule
28   	 * @author Mark Turner
29   	 *
30   	 * The approach is based on the following papers.
31   	 *
32   	 * M. Turner, T. Berthold, M. Besancon, T. Koch@n
33   	 * Branching via Cutting Plane Selection: Improving Hybrid Branching,@n
34   	 * arXiv preprint arXiv:2306.06050
35   	 *
36   	 * The Gomory cut branching rule selects a candidate integer variable $j$ with a fractional solution value.
37   	 * Each candidate variable must be a basic variable in the LP Tableau (if not then it would have to be at its bound
38   	 * that is integer-valued)
39   	 * This branching rule calculates the GMI cut for the aggregated row of the LP tableau associated with the
40   	 * candidate variable.
41   	 * The generated cut is then scored using a weighted sum rule.
42   	 * The branching candidate whose cut is highest scoring is then selected.
43   	 * For more details on the method, see:
44   	 *
45   	 * @par
46   	 * Mark Turner, Timo Berthold, Mathieu Besançon, Thorsten Koch@n
47   	 * Branching via Cutting Plane Selection: Improving Hybrid Branching@n
48   	 * 2023@n
49   	 */
50   	
51   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
52   	
53   	#include "scip/branch_gomory.h"
54   	#include "scip/pub_branch.h"
55   	#include "scip/pub_var.h"
56   	#include "scip/pub_lp.h"
57   	#include "scip/pub_tree.h"
58   	#include "scip/scip_branch.h"
59   	#include "scip/scip_cut.h"
60   	#include "scip/scip_mem.h"
61   	#include "scip/scip_message.h"
62   	#include "scip/scip_numerics.h"
63   	#include "scip/scip_lp.h"
64   	#include "scip/scip_tree.h"
65   	#include "scip_param.h"
66   	#include "scip/branch_relpscost.h"
67   	#include <string.h>
68   	#include <assert.h>
69   	
70   	
71   	
72   	#define BRANCHRULE_NAME            "gomory"
73   	#define BRANCHRULE_DESC            "Gomory cut score branching"
74   	#define BRANCHRULE_PRIORITY        -1000
75   	#define BRANCHRULE_MAXDEPTH        -1
76   	#define BRANCHRULE_MAXBOUNDDIST    1.0
77   	
78   	#define DEFAULT_MAXNCANDS          -1    /**< maximum number of branching candidates to produce a cut for */
79   	#define DEFAULT_EFFICACYWEIGHT     1.0   /**< the weight of efficacy in weighted sum cut scoring rule */
80   	#define DEFAULT_OBJPARALLELWEIGHT  0.0   /**< the weight of objective parallelism in weighted sum scoring rule */
81   	#define DEFAULT_INTSUPPORTWEIGHT   0.0   /**< the weight of integer support in weighted sum cut scoring rule */
82   	#define DEFAULT_PERFORMRELPSCOST   FALSE /**< if relpscost branching should be called without actual branching */
83   	#define DEFAULT_USEWEAKERCUTS      TRUE  /**< use weaker cuts derived from the exact branching split */
84   	
85   	
86   	/*
87   	 * Data structures
88   	 */
89   	
90   	/** branching rule data */
91   	struct SCIP_BranchruleData
92   	{
93   	   int                   maxncands;             /**< maximum number of variable candidates to produce cut for */
94   	   SCIP_Real             efficacyweight;        /**< the weight of efficacy in weighted sum cut scoring rule */
95   	   SCIP_Real             objparallelweight;     /**< the weight of objective parallelism in weighted sum scoring rule */
96   	   SCIP_Real             intsupportweight;      /**< the weight of integer support in weighted sum cut scoring rule */
97   	   SCIP_Bool             performrelpscost;      /**< if relpscost branching should be called without actual branching */
98   	   SCIP_Bool             useweakercuts;         /**< use weaker cuts derived from the exact branching split */
99   	};
100  	
101  	
102  	/*
103  	 * Local methods
104  	 */
105  	
106  	/** Generate GMI cut: The GMI is given by
107  	    * sum(f_j x_j                  , j in J_I s.t. f_j <= f_0) +
108  	    * sum((1-f_j)*f_0/(1 - f_0) x_j, j in J_I s.t. f_j  > f_0) +
109  	    * sum(a_j x_j,                 , j in J_C s.t. a_j >=   0) -
110  	    * sum(a_j*f_0/(1-f_0) x_j      , j in J_C s.t. a_j  <   0) >= f_0.
111  	    * where J_I are the integer non-basic variables and J_C are the continuous.
112  	    * f_0 is the fractional part of lpval
113  	    * a_j is the j-th coefficient of the tableau row and f_j its fractional part
114  	    * Note: we create -% <= -f_0 !!
115  	    * Note: this formula is valid for a problem of the form Ax = b, x>= 0. Since we do not have
116  	    * such problem structure in general, we have to (implicitly) transform whatever we are given
117  	    * to that form. Specifically, non-basic variables at their lower bound are shifted so that the lower
118  	    * bound is 0 and non-basic at their upper bound are complemented. */
119  	static
120  	SCIP_Bool getGMIFromRow(
121  	   SCIP*                 scip,               /**< SCIP data structure */
122  	   int                   ncols,              /**< Number of columns (original variables) in the LP */
123  	   int                   nrows,              /**< Number of rows (slack variables) in the LP */
124  	   SCIP_COL**            cols,               /**< Column data of the LP */
125  	   SCIP_ROW**            rows,               /**< Row data of the LP */
126  	   const SCIP_Real*      binvrow,            /**< row of B^-1 for current basic variable */
127  	   const SCIP_Real*      binvarow,           /**< row of B^-1A for current basic variable */
128  	   const SCIP_Real*      lpval,              /**< value of variable at current LP solution */
129  	   SCIP_Real*            cutcoefs,           /**< array to store cut coefficients */
130  	   SCIP_Real*            cutrhs,             /**< pointer to store rhs of cut */
131  	   SCIP_Bool             useweakerscuts      /**< use weakener cuts derived from the exact branching split */
132  	   )
133  	{
134  	   SCIP_COL** rowcols;
135  	   SCIP_COL* col;
136  	   SCIP_ROW* row;
137  	   SCIP_Real* rowvals;
138  	   SCIP_BASESTAT basestat;
139  	   SCIP_Real rowelem;
140  	   SCIP_Real cutelem;
141  	   SCIP_Real f0;
142  	   SCIP_Real f0complementratio;
143  	   SCIP_Real rowrhs;
144  	   SCIP_Real rowlhs;
145  	   SCIP_Real rowact;
146  	   SCIP_Real rowrhsslack;
147  	   int i;
148  	   int c;
149  	
150  	   /* Clear the memory array of cut coefficients. It may store that of the last computed cut */
151  	   BMSclearMemoryArray(cutcoefs, ncols);
152  	
153  	   /* compute fractionality f0 and f0/(1-f0) */
154  	   f0 = SCIPfeasFrac(scip, *lpval);
155  	   f0complementratio = f0 / (1.0 - f0);
156  	
157  	   /* The rhs of the cut is the fractional part of the LP solution of the basic variable */
158  	   *cutrhs = -f0;
159  	
160  	   /* Generate cut coefficient for the original variables */
161  	   for ( c = 0; c < ncols; c++ )
162  	   {
163  	      col = cols[c];
164  	      assert( col != NULL );
165  	
166  	      basestat = SCIPcolGetBasisStatus(col);
167  	      /* Get simplex tableau coefficient */
168  	      if ( basestat == SCIP_BASESTAT_LOWER )
169  	      {
170  	         /* Take coefficient if nonbasic at lower bound */
171  	         rowelem = binvarow[c];
172  	      }
173  	      else if ( basestat == SCIP_BASESTAT_UPPER )
174  	      {
175  	         /* Flip coefficient if nonbasic at upper bound: x --> u - x */
176  	         rowelem = -binvarow[c];
177  	      }
178  	      else
179  	      {
180  	         /* Nonbasic free variable at zero or basic variable. Just skip it. */
181  	         continue;
182  	      }
183  	
184  	      /* Integer variables */
185  	      if ( SCIPcolIsIntegral(col) && !useweakerscuts )
186  	      {
187  	         /* Warning: Because of numerics we can have cutelem < 0
188  	          * In such a case it is very close to 0, so isZero will catch and we can ignore the coefficient */
189  	         cutelem = SCIPfrac(scip, rowelem);
190  	         if ( cutelem > f0 )
191  	         {
192  	            /* sum((1 - f_j) * f_0/(1 - f_0) x_j, j in J_I s.t. f_j > f_0 */
193  	            cutelem = -((1.0 - cutelem) * f0complementratio);
194  	         }
195  	         else
196  	         {
197  	            /* sum(f_j * x_j, j in J_I s.t. f_j <= 0 */
198  	            cutelem = -cutelem;
199  	         }
200  	      }
201  	      /* Then continuous variables */
202  	      else
203  	      {
204  	         if ( rowelem < 0 )
205  	         {
206  	            /* sum(a_j* f_0/(1 - f_0) x_j, j in J_C s.t. a_j < 0 */
207  	            cutelem = rowelem * f0complementratio;
208  	         }
209  	         else
210  	         {
211  	            /* sum(a_j * x_j, j in J_C s.t. a_j >= 0 */
212  	            cutelem = -rowelem;
213  	         }
214  	      }
215  	
216  	      /* Cut is defined when variables are in [0, infinity). Translate to general bounds. */
217  	      if ( !SCIPisZero(scip, cutelem) )
218  	      {
219  	         if ( basestat == SCIP_BASESTAT_UPPER )
220  	         {
221  	            cutelem = -cutelem;
222  	            *cutrhs += cutelem * SCIPcolGetUb(col);
223  	         }
224  	         else
225  	         {
226  	            *cutrhs += cutelem * SCIPcolGetLb(col);
227  	         }
228  	         /* Add coefficient to cut */
229  	         cutcoefs[SCIPcolGetLPPos(col)] = cutelem;
230  	      }
231  	   }
232  	
233  	   /* generate cut coefficient for the slack variables. Skip the basic ones */
234  	   for ( c = 0; c < nrows; c++ )
235  	   {
236  	      row = rows[c];
237  	      assert( row != NULL );
238  	      basestat = SCIProwGetBasisStatus(row);
239  	
240  	      /* Get the simplex tableau coefficient */
241  	      if ( basestat == SCIP_BASESTAT_LOWER )
242  	      {
243  	         /* Take coefficient if nonbasic at lower bound */
244  	         rowelem = binvrow[SCIProwGetLPPos(row)];
245  	         /* If there is a >= constraint or ranged constraint at the lower bound, we have to flip the row element */
246  	         if ( !SCIPisInfinity(scip, -SCIProwGetLhs(row)) )
247  	            rowelem = -rowelem;
248  	      }
249  	      else if ( basestat == SCIP_BASESTAT_UPPER )
250  	      {
251  	         /* Take element if nonbasic at upper bound. Only non-positive slack variables can be nonbasic at upper,
252  	          * therefore they should be flipped twice, meaning we can take the element directly */
253  	         rowelem = binvrow[SCIProwGetLPPos(row)];
254  	      }
255  	      else
256  	      {
257  	         /* Nonbasic free variable at zero or basic variable. Free variable should not happen here. Just skip if free */
258  	         assert( basestat == SCIP_BASESTAT_BASIC );
259  	         continue;
260  	      }
261  	
262  	      /* Integer rows */
263  	      if ( SCIProwIsIntegral(row) && !SCIProwIsModifiable(row) && !useweakerscuts )
264  	      {
265  	         /* Warning: Because of numerics we can have cutelem < 0
266  	         * In such a case it is very close to 0, so isZero will catch and we can ignore the coefficient */
267  	         cutelem = SCIPfrac(scip, rowelem);
268  	         if ( cutelem > f0 )
269  	         {
270  	            /* sum((1 - f_j) * f_0/(1 - f_0) x_j, j in J_I s.t. f_j > f_0 */
271  	            cutelem = -((1.0 - cutelem) * f0complementratio);
272  	         }
273  	         else
274  	         {
275  	            /* sum(f_j * x_j, j in J_I s.t. f_j <= 0 */
276  	            cutelem = -cutelem;
277  	         }
278  	      }
279  	      /* Then continuous variables */
280  	      else
281  	      {
282  	         if ( rowelem < 0 )
283  	         {
284  	            /* sum(a_j* f_0/(1 - f_0) x_j, j in J_C s.t. a_j < 0 */
285  	            cutelem = rowelem * f0complementratio;
286  	         }
287  	         else
288  	         {
289  	            /* sum(a_j * x_j, j in J_C s.t. a_j >= 0 */
290  	            cutelem = -rowelem;
291  	         }
292  	      }
293  	
294  	      /* Cut is defined in original variables, so we replace slack variables by their original definition */
295  	      if ( !SCIPisZero(scip, cutelem) )
296  	      {
297  	         /* Coefficient is large enough so we keep it */
298  	         rowlhs = SCIProwGetLhs(row);
299  	         rowrhs = SCIProwGetRhs(row);
300  	         assert( SCIPisLE(scip, rowlhs, rowrhs) );
301  	         assert( !SCIPisInfinity(scip, rowlhs) || !SCIPisInfinity(scip, rowrhs) );
302  	
303  	         /* If the slack variable is fixed we can ignore this cut coefficient */
304  	         if ( SCIPisFeasZero(scip, rowrhs - rowlhs) )
305  	            continue;
306  	
307  	         /* Un-flip sack variable and adjust rhs if necessary.
308  	          * Row at lower basis means the slack variable is at its upper bound.
309  	          * Since SCIP adds +1 slacks, this can only happen when constraints have finite lhs */
310  	         if ( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER )
311  	         {
312  	            assert( !SCIPisInfinity(scip, -rowlhs) );
313  	            cutelem = -cutelem;
314  	         }
315  	
316  	         rowcols = SCIProwGetCols(row);
317  	         rowvals = SCIProwGetVals(row);
318  	
319  	         /* Eliminate slack variables. rowcols is sorted [columns in LP, columns not in LP] */
320  	         for ( i = 0; i < SCIProwGetNLPNonz(row); i++ )
321  	            cutcoefs[SCIPcolGetLPPos(rowcols[i])] -= cutelem * rowvals[i];
322  	
323  	         /* Modify the rhs */
324  	         rowact = SCIPgetRowActivity(scip, row);
325  	         rowrhsslack = rowrhs - rowact;
326  	
327  	         if ( SCIPisFeasZero(scip, rowrhsslack) )
328  	            *cutrhs -= cutelem * (rowrhs - SCIProwGetConstant(row));
329  	         else
330  	         {
331  	            assert( SCIPisFeasZero(scip, rowact - rowlhs) );
332  	            *cutrhs -= cutelem * (rowlhs - SCIProwGetConstant(row));
333  	         }
334  	
335  	      }
336  	   }
337  	
338  	   return TRUE;
339  	}
340  	
341  	
342  	/*
343  	 * Callback methods of branching rule
344  	 */
345  	
346  	
347  	/** copy method for branchrule plugins (called when SCIP copies plugins) */
348  	static
349  	SCIP_DECL_BRANCHCOPY(branchCopyGomory)
350  	{  /*lint --e{715}*/
351  	   assert(scip != NULL);
352  	   assert(branchrule != NULL);
353  	   assert(strcmp(SCIPbranchruleGetName(branchrule), BRANCHRULE_NAME) == 0);
354  	
355  	   /* call inclusion method of branchrule */
356  	   SCIP_CALL( SCIPincludeBranchruleGomory(scip) );
357  	
358  	   return SCIP_OKAY;
359  	}
360  	
361  	
362  	/** destructor of branching rule to free user data (called when SCIP is exiting) */
363  	static
364  	SCIP_DECL_BRANCHFREE(branchFreeGomory)
365  	{  /*lint --e{715}*/
366  	   SCIP_BRANCHRULEDATA* branchruledata;
367  	
368  	   /* free branching rule data */
369  	   branchruledata = SCIPbranchruleGetData(branchrule);
370  	   assert(branchruledata != NULL);
371  	
372  	   SCIPfreeBlockMemoryNull(scip, &branchruledata);
373  	
374  	   return SCIP_OKAY;
375  	}
376  	
377  	
378  	/** branching execution method for fractional LP solutions */
379  	static
380  	SCIP_DECL_BRANCHEXECLP(branchExeclpGomory)
381  	{  /*lint --e{715}*/
382  	   SCIP_BRANCHRULEDATA* branchruledata;
383  	   SCIP_VAR** lpcands;
384  	   SCIP_COL** cols;
385  	   SCIP_ROW** rows;
386  	   SCIP_Real* lpcandssol;
387  	   SCIP_Real* lpcandsfrac;
388  	   SCIP_Real* binvrow;
389  	   SCIP_Real* binvarow;
390  	   SCIP_Real* cutcoefs;
391  	   SCIP_ROW* cut;
392  	   SCIP_COL* col;
393  	   int* basisind;
394  	   int* basicvarpos2tableaurow;
395  	   int* inds;
396  	   const char* name;
397  	   SCIP_Real cutrhs;
398  	   SCIP_Real score;
399  	   SCIP_Real bestscore;
400  	   SCIP_Bool success;
401  	   int nlpcands;
402  	   int maxncands;
403  	   int ncols;
404  	   int nrows;
405  	   int lppos;
406  	   int ninds;
407  	   int bestcand;
408  	
409  	   name = (char *) "test";
410  	
411  	   assert(branchrule != NULL);
412  	   assert(strcmp(SCIPbranchruleGetName(branchrule), BRANCHRULE_NAME) == 0);
413  	   assert(scip != NULL);
414  	   assert(result != NULL);
415  	
416  	   SCIPdebugMsg(scip, "Execlp method of Gomory branching in node %" SCIP_LONGINT_FORMAT "\n", SCIPnodeGetNumber(SCIPgetCurrentNode(scip)));
417  	
418  	   if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
419  	   {
420  	      *result = SCIP_DIDNOTRUN;
421  	      SCIPdebugMsg(scip, "Could not apply Gomory branching, as the current LP was not solved to optimality.\n");
422  	
423  	      return SCIP_OKAY;
424  	   }
425  	
426  	   /* Get branching candidates */
427  	   SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, &lpcandsfrac, NULL, &nlpcands, NULL) );
428  	   assert(nlpcands > 0);
429  	
430  	   *result = SCIP_DIDNOTRUN;
431  	
432  	   /* Get branching rule data */
433  	   branchruledata = SCIPbranchruleGetData(branchrule);
434  	   assert(branchruledata != NULL);
435  	
436  	   /* Compute the reliability pseudo-cost branching scores for the candidates */
437  	   if ( branchruledata->performrelpscost )
438  	   {
439  	      /* We do not branch using this rule, but if enabled do take all the bound and conflict inferences made */
440  	      SCIP_CALL( SCIPexecRelpscostBranching(scip, lpcands, lpcandssol, lpcandsfrac, nlpcands, FALSE,  result) );
441  	      assert(*result == SCIP_DIDNOTRUN || *result == SCIP_CUTOFF || *result == SCIP_REDUCEDDOM);
442  	   }
443  	
444  	   /* Return SCIP_OKAY if relpscost has shown that this node can be cutoff or some variable domains have changed */
445  	   if( *result == SCIP_CUTOFF || *result == SCIP_REDUCEDDOM )
446  	   {
447  	      return SCIP_OKAY;
448  	   }
449  	
450  	   /* Get the maximum number of LP branching candidates that we generate cuts for and score */
451  	   if( branchruledata->maxncands >= 0 )
452  	   {
453  	      maxncands = MIN(nlpcands, branchruledata->maxncands);
454  	   }
455  	   else
456  	   {
457  	      maxncands = nlpcands;
458  	   }
459  	
460  	   /* Get the Column and Row data */
461  	   SCIP_CALL(SCIPgetLPColsData(scip, &cols, &ncols));
462  	   SCIP_CALL(SCIPgetLPRowsData(scip, &rows, &nrows));
463  	
464  	   /* Allocate temporary memory */
465  	   SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, ncols) );
466  	   SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) );
467  	   SCIP_CALL( SCIPallocBufferArray(scip, &basicvarpos2tableaurow, ncols) );
468  	   SCIP_CALL( SCIPallocBufferArray(scip, &binvrow, nrows) );
469  	   SCIP_CALL( SCIPallocBufferArray(scip, &binvarow, ncols) );
470  	   SCIP_CALL( SCIPallocBufferArray(scip, &inds, nrows) );
471  	
472  	   /* Create basis indices mapping (from the column position to LP tableau rox index) */
473  	   for( int i = 0; i < ncols; ++i )
474  	   {
475  	      basicvarpos2tableaurow[i] = -1;
476  	   }
477  	   SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) );
478  	   for( int i = 0; i < nrows; ++i )
479  	   {
480  	      if( basisind[i] >= 0 )
481  	         basicvarpos2tableaurow[basisind[i]] = i;
482  	   }
483  	
484  	   /* Initialise the best candidate */
485  	   bestcand = 0;
486  	   bestscore = -SCIPinfinity(scip);
487  	   ninds = -1;
488  	
489  	   /* Iterate over candidates and get best cut score */
490  	   for( int i = 0; i < maxncands; i++ ) {
491  	
492  	      /* Initialise the score of the cut */
493  	      score = 0;
494  	
495  	      /* Get the LP position of the branching candidate */
496  	      col = SCIPvarGetCol(lpcands[i]);
497  	      lppos = SCIPcolGetLPPos(col);
498  	      assert(lppos != -1);
499  	
500  	      /* get the row of B^-1 for this basic integer variable with fractional solution value */
501  	      SCIP_CALL(SCIPgetLPBInvRow(scip, basicvarpos2tableaurow[lppos], binvrow, inds, &ninds));
502  	
503  	      /* Get the Tableau row for this basic integer variable with fractional solution value */
504  	      SCIP_CALL(SCIPgetLPBInvARow(scip, basicvarpos2tableaurow[lppos], binvrow, binvarow, inds, &ninds));
505  	
506  	      /* Compute the GMI cut */
507  	      success = getGMIFromRow(scip, ncols, nrows, cols, rows, binvrow, binvarow, &lpcandssol[i], cutcoefs,
508  	         &cutrhs, branchruledata->useweakercuts);
509  	
510  	      /* Calculate the weighted sum score of measures */
511  	      if ( success )
512  	      {
513  	         cut = NULL;
514  	         SCIP_CALL( SCIPcreateEmptyRowUnspec(scip, &cut, name, -SCIPinfinity(scip), cutrhs, TRUE,
515  	                   FALSE, TRUE) );
516  	         for( int j = 0; j < ncols; ++j )
517  	         {
518  	            if( !SCIPisZero(scip, cutcoefs[j]) )
519  	            {
520  	               SCIP_CALL( SCIPaddVarToRow(scip, cut, SCIPcolGetVar(cols[j]),
521  	                                          cutcoefs[SCIPcolGetLPPos(cols[j])]) );
522  	            }
523  	         }
524  	         assert( SCIPgetCutEfficacy(scip, NULL, cut) >= -SCIPfeastol(scip) );
525  	         if ( branchruledata-> efficacyweight != 0 )
526  	            score += branchruledata->efficacyweight * SCIPgetCutEfficacy(scip, NULL, cut);
527  	         if ( branchruledata->objparallelweight != 0 )
528  	            score += branchruledata->objparallelweight * SCIPgetRowObjParallelism(scip, cut);
529  	         if ( branchruledata->intsupportweight != 0 )
530  	            score += branchruledata->intsupportweight * SCIPgetRowNumIntCols(scip, cut) / (SCIP_Real) SCIProwGetNNonz(cut);
531  	         SCIP_CALL( SCIPreleaseRow(scip, &cut) );
532  	
533  	         /* Replace the best cut if score is higher */
534  	         if (score > bestscore) {
535  	            bestscore = score;
536  	            bestcand = i;
537  	         }
538  	      }
539  	   }
540  	
541  	   /* Free temporary memory */
542  	   SCIPfreeBufferArray(scip, &inds);
543  	   SCIPfreeBufferArray(scip, &binvrow);
544  	   SCIPfreeBufferArray(scip, &binvarow);
545  	   SCIPfreeBufferArray(scip, &basicvarpos2tableaurow);
546  	   SCIPfreeBufferArray(scip, &basisind);
547  	   SCIPfreeBufferArray(scip, &cutcoefs);
548  	
549  	   SCIPdebugMsg(scip, " -> %d candidates, selected candidate %d: variable <%s> (frac=%g, factor=%g, score=%g)\n",
550  	                nlpcands, bestcand, SCIPvarGetName(lpcands[bestcand]), lpcandsfrac[bestcand],
551  	                SCIPvarGetBranchFactor(lpcands[bestcand]), bestscore);
552  	
553  	   /* Perform the branching */
554  	   SCIP_CALL( SCIPbranchVar(scip, lpcands[bestcand], NULL, NULL, NULL) );
555  	   *result = SCIP_BRANCHED;
556  	
557  	   return SCIP_OKAY;
558  	}
559  	
560  	
561  	/*
562  	 * branching rule specific interface methods
563  	 */
564  	
565  	/** creates the Gomory cut branching rule and includes it in SCIP */
566  	SCIP_RETCODE SCIPincludeBranchruleGomory(
567  	   SCIP*                 scip                /**< SCIP data structure */
568  	)
569  	{
570  	   SCIP_BRANCHRULEDATA* branchruledata;
571  	   SCIP_BRANCHRULE* branchrule;
572  	
573  	   /* create branching rule data */
574  	   SCIP_CALL( SCIPallocBlockMemory(scip, &branchruledata) );
575  	
576  	   /* include branching rule */
577  	   SCIP_CALL( SCIPincludeBranchruleBasic(scip, &branchrule, BRANCHRULE_NAME, BRANCHRULE_DESC, BRANCHRULE_PRIORITY,
578  	                                         BRANCHRULE_MAXDEPTH, BRANCHRULE_MAXBOUNDDIST, branchruledata) );
579  	
580  	   assert(branchrule != NULL);
581  	
582  	   /* set non-fundamental callbacks via specific setter functions*/
583  	   SCIP_CALL( SCIPsetBranchruleCopy(scip, branchrule, branchCopyGomory) );
584  	   SCIP_CALL( SCIPsetBranchruleFree(scip, branchrule, branchFreeGomory) );
585  	   SCIP_CALL( SCIPsetBranchruleExecLp(scip, branchrule, branchExeclpGomory) );
586  	
587  	   /* Gomory cut branching rule parameters */
588  	   SCIP_CALL( SCIPaddIntParam(scip,"branching/gomory/maxncands",
589  	         "maximum amount of branching candidates to generate Gomory cuts for (-1: all candidates)",
590  	         &branchruledata->maxncands, FALSE, DEFAULT_MAXNCANDS, -1, INT_MAX, NULL, NULL) );
591  	   SCIP_CALL( SCIPaddRealParam(scip,"branching/gomory/efficacyweight",
592  	         "weight of efficacy in the weighted sum cut scoring rule",
593  	         &branchruledata->efficacyweight, FALSE, DEFAULT_EFFICACYWEIGHT, -1.0, 1.0, NULL, NULL) );
594  	   SCIP_CALL( SCIPaddRealParam(scip,"branching/gomory/objparallelweight",
595  	         "weight of objective parallelism in the weighted sum cut scoring rule",
596  	         &branchruledata->objparallelweight, FALSE, DEFAULT_OBJPARALLELWEIGHT, -1.0, 1.0, NULL, NULL) );
597  	   SCIP_CALL( SCIPaddRealParam(scip,"branching/gomory/intsupportweight",
598  	         "weight of integer support in the weighted sum cut scoring rule",
599  	         &branchruledata->intsupportweight, FALSE, DEFAULT_INTSUPPORTWEIGHT, -1.0, 1.0, NULL, NULL) );
600  	   SCIP_CALL( SCIPaddBoolParam(scip,"branching/gomory/performrelpscost",
601  	         "whether relpscost branching should be called without branching (used for bound inferences and conflicts)",
602  	         &branchruledata->performrelpscost, FALSE, DEFAULT_PERFORMRELPSCOST, NULL, NULL) );
603  	   SCIP_CALL( SCIPaddBoolParam(scip,"branching/gomory/useweakercuts",
604  	         "use weaker cuts that are exactly derived from the branching split disjunction",
605  	         &branchruledata->useweakercuts, FALSE, DEFAULT_USEWEAKERCUTS, NULL, NULL) );
606  	
607  	   return SCIP_OKAY;
608  	}
609