1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the program and library             */
4    	/*         SCIP --- Solving Constraint Integer Programs                      */
5    	/*                                                                           */
6    	/*  Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB)                      */
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   sepa_gomory.c
26   	 * @ingroup DEFPLUGINS_SEPA
27   	 * @brief  Gomory MIR Cuts
28   	 * @author Tobias Achterberg
29   	 * @author Stefan Heinz
30   	 * @author Domenico Salvagnin
31   	 * @author Marc Pfetsch
32   	 * @author Leona Gottwald
33   	 */
34   	
35   	/**@todo try k-Gomory-cuts (s. Cornuejols: K-Cuts: A Variation of Gomory Mixed Integer Cuts from the LP Tableau)
36   	 *
37   	 * @todo Try cuts on the objective tableau row.
38   	 *
39   	 * @todo Also try negative basis inverse row?
40   	 *
41   	 * @todo It happens that the SCIPcalcMIR() function returns with the same cut for different calls. Check if this is a
42   	 *       bug or do not use it for the MIP below and turn off presolving and all heuristics:
43   	 *
44   	 *  Max y
45   	 *  Subject to
46   	 *  c1: -x + y <= 1
47   	 *  c2: 2x + 3y <= 12
48   	 *  c3: 3x + 2y <= 12
49   	 *  Bounds
50   	 *  0 <= x
51   	 *  0 <= y
52   	 *  General
53   	 *  x
54   	 *  y
55   	 *  END
56   	 */
57   	
58   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
59   	
60   	#include "blockmemshell/memory.h"
61   	#include "scip/cuts.h"
62   	#include "scip/pub_lp.h"
63   	#include "scip/pub_message.h"
64   	#include "scip/pub_misc.h"
65   	#include "scip/pub_misc_sort.h"
66   	#include "scip/pub_sepa.h"
67   	#include "scip/pub_var.h"
68   	#include "scip/scip_branch.h"
69   	#include "scip/scip_cut.h"
70   	#include "scip/scip_general.h"
71   	#include "scip/scip_lp.h"
72   	#include "scip/scip_mem.h"
73   	#include "scip/scip_message.h"
74   	#include "scip/scip_numerics.h"
75   	#include "scip/scip_param.h"
76   	#include "scip/scip_prob.h"
77   	#include "scip/scip_randnumgen.h"
78   	#include "scip/scip_sepa.h"
79   	#include "scip/scip_solvingstats.h"
80   	#include "scip/scip_tree.h"
81   	#include "scip/scip_var.h"
82   	#include "scip/sepa_gomory.h"
83   	#include "struct_history.h"
84   	#include <string.h>
85   	
86   	#define SEPA_NAME              "gomory"
87   	#define SEPA_DESC              "separator for Gomory mixed-integer and strong CG cuts from LP tableau rows"
88   	#define SEPA_PRIORITY             -1000
89   	#define SEPA_FREQ                    10
90   	#define SEPA_MAXBOUNDDIST           1.0
91   	#define SEPA_USESSUBSCIP          FALSE /**< does the separator use a secondary SCIP instance? */
92   	#define SEPA_DELAY                FALSE /**< should separation method be delayed, if other separators found cuts? */
93   	
94   	#define DEFAULT_MAXROUNDS             5 /**< maximal number of gomory separation rounds per node (-1: unlimited) */
95   	#define DEFAULT_MAXROUNDSROOT        10 /**< maximal number of gomory separation rounds in the root node (-1: unlimited) */
96   	#define DEFAULT_MAXSEPACUTS          50 /**< maximal number of gomory cuts separated per separation round */
97   	#define DEFAULT_MAXSEPACUTSROOT     200 /**< maximal number of gomory cuts separated per separation round in root node */
98   	#define DEFAULT_MAXRANK              -1 /**< maximal rank of a gomory cut that could not be scaled to integral coefficients (-1: unlimited) */
99   	#define DEFAULT_MAXRANKINTEGRAL      -1 /**< maximal rank of a gomory cut that could be scaled to integral coefficients (-1: unlimited) */
100  	#define DEFAULT_DYNAMICCUTS        TRUE /**< should generated cuts be removed from the LP if they are no longer tight? */
101  	#define DEFAULT_AWAY               0.01 /**< minimal integrality violation of a basis variable in order to try Gomory cut */
102  	#define DEFAULT_MAKEINTEGRAL      FALSE /**< try to scale all cuts to integral coefficients */
103  	#define DEFAULT_FORCECUTS          TRUE /**< if conversion to integral coefficients failed still consider the cut */
104  	#define DEFAULT_SEPARATEROWS       TRUE /**< separate rows with integral slack */
105  	#define DEFAULT_DELAYEDCUTS       FALSE /**< should cuts be added to the delayed cut pool? */
106  	#define DEFAULT_SIDETYPEBASIS      TRUE /**< choose side types of row (lhs/rhs) based on basis information? */
107  	#define DEFAULT_TRYSTRONGCG        TRUE /**< try to generate strengthened Chvatal-Gomory cuts? */
108  	#define DEFAULT_GENBOTHGOMSCG      TRUE /**< should both Gomory and strong CG cuts be generated (otherwise take best) */
109  	#define DEFAULT_RANDSEED             53 /**< initial random seed */
110  	
111  	#define BOUNDSWITCH              0.9999 /**< threshold for bound switching - see SCIPcalcMIR() */
112  	#define POSTPROCESS                TRUE /**< apply postprocessing after MIR calculation - see SCIPcalcMIR() */
113  	#define USEVBDS                    TRUE /**< use variable bounds - see SCIPcalcMIR() */
114  	#define FIXINTEGRALRHS            FALSE /**< try to generate an integral rhs - see SCIPcalcMIR() */
115  	#define MAKECONTINTEGRAL          FALSE /**< convert continuous variable to integral variables in SCIPmakeRowIntegral() */
116  	
117  	#define MAXAGGRLEN(nvars)          (0.1*(nvars)+1000) /**< maximal length of base inequality */
118  	
119  	
120  	/** separator data */
121  	struct SCIP_SepaData
122  	{
123  	   SCIP_RANDNUMGEN*      randnumgen;         /**< random number generator */
124  	   SCIP_SEPA*            strongcg;           /**< strong CG cut separator */
125  	   SCIP_SEPA*            gomory;             /**< gomory cut separator */
126  	   SCIP_Real             away;               /**< minimal integrality violation of a basis variable in order to try Gomory cut */
127  	   int                   maxrounds;          /**< maximal number of gomory separation rounds per node (-1: unlimited) */
128  	   int                   maxroundsroot;      /**< maximal number of gomory separation rounds in the root node (-1: unlimited) */
129  	   int                   maxsepacuts;        /**< maximal number of gomory cuts separated per separation round */
130  	   int                   maxsepacutsroot;    /**< maximal number of gomory cuts separated per separation round in root node */
131  	   int                   maxrank;            /**< maximal rank of a gomory cut that could not be scaled to integral coefficients (-1: unlimited) */
132  	   int                   maxrankintegral;    /**< maximal rank of a gomory cut that could be scaled to integral coefficients (-1: unlimited) */
133  	   int                   lastncutsfound;     /**< total number of cuts found after last call of separator */
134  	   SCIP_Bool             dynamiccuts;        /**< should generated cuts be removed from the LP if they are no longer tight? */
135  	   SCIP_Bool             makeintegral;       /**< try to scale all cuts to integral coefficients */
136  	   SCIP_Bool             forcecuts;          /**< if conversion to integral coefficients failed still consider the cut */
137  	   SCIP_Bool             separaterows;       /**< separate rows with integral slack */
138  	   SCIP_Bool             delayedcuts;        /**< should cuts be added to the delayed cut pool? */
139  	   SCIP_Bool             sidetypebasis;      /**< choose side types of row (lhs/rhs) based on basis information? */
140  	   SCIP_Bool             trystrongcg;        /**< try to generate strengthened Chvatal-Gomory cuts? */
141  	   SCIP_Bool             genbothgomscg;      /**< should both Gomory and strong CG cuts be generated (otherwise take best) */
142  	};
143  	
144  	
145  	/** returns TRUE if the cut can be taken, otherwise FALSE if there some numerical evidences */
146  	static
147  	SCIP_RETCODE evaluateCutNumerics(
148  	   SCIP*                 scip,               /**< SCIP data structure */
149  	   SCIP_SEPADATA*        sepadata,           /**< data of the separator */
150  	   SCIP_ROW*             cut,                /**< cut to check */
151  	   SCIP_Longint          maxdnom,            /**< maximal denominator to use for scaling */
152  	   SCIP_Real             maxscale,           /**< maximal scaling factor */
153  	   SCIP_Bool*            useful              /**< pointer to store if the cut is useful */
154  	   )
155  	{
156  	   SCIP_Bool madeintegral = FALSE;
157  	
158  	   assert(useful != NULL);
159  	
160  	   *useful = FALSE;
161  	
162  	   if( sepadata->makeintegral && SCIPgetRowNumIntCols(scip, cut) == SCIProwGetNNonz(cut) )
163  	   {
164  	      /* try to scale the cut to integral values */
165  	      SCIP_CALL( SCIPmakeRowIntegral(scip, cut, -SCIPepsilon(scip), SCIPsumepsilon(scip),
166  	            maxdnom, maxscale, MAKECONTINTEGRAL, &madeintegral) );
167  	
168  	      if( !madeintegral && !sepadata->forcecuts )
169  	         return SCIP_OKAY;
170  	
171  	      /* in case the right hand side is plus infinity (due to scaling) the cut is useless so we are not taking it at all */
172  	      if( madeintegral && SCIPisInfinity(scip, SCIProwGetRhs(cut)) )
173  	         return SCIP_OKAY;
174  	   }
175  	
176  	   /* discard integral cut if the rank is too high */
177  	   if( madeintegral && sepadata->maxrankintegral != -1 && (SCIProwGetRank(cut) > sepadata->maxrankintegral) )
178  	      return SCIP_OKAY;
179  	
180  	   /* discard cut if the rank is too high */
181  	   if( !madeintegral && (sepadata->maxrank != -1) && (SCIProwGetRank(cut) > sepadata->maxrank) )
182  	      return SCIP_OKAY;
183  	
184  	   *useful = TRUE;
185  	
186  	   return SCIP_OKAY;
187  	}
188  	
189  	
190  	/** add cut */
191  	static
192  	SCIP_RETCODE addCut(
193  	   SCIP*                 scip,               /**< SCIP instance */
194  	   SCIP_SEPADATA*        sepadata,           /**< separator data */
195  	   SCIP_VAR**            vars,               /**< array of variables */
196  	   int                   c,                  /**< index of basic variable (< 0 for slack variables) */
197  	   SCIP_Longint          maxdnom,            /**< maximal denominator to use for scaling */
198  	   SCIP_Real             maxscale,           /**< maximal scaling factor */
199  	   int                   cutnnz,             /**< number of nonzeros in cut */
200  	   int*                  cutinds,            /**< variable indices in cut */
201  	   SCIP_Real*            cutcoefs,           /**< cut cofficients */
202  	   SCIP_Real             cutefficacy,        /**< cut efficacy */
203  	   SCIP_Real             cutrhs,             /**< rhs of cut */
204  	   SCIP_Bool             cutislocal,         /**< whether cut is local */
205  	   int                   cutrank,            /**< rank of cut */
206  	   SCIP_Bool             strongcg,           /**< whether the cut arises from the strong-CG procedure */
207  	   SCIP_Bool*            cutoff,             /**< pointer to store whether a cutoff appeared */
208  	   int*                  naddedcuts          /**< pointer to store number of added cuts */
209  	   )
210  	{
211  	   assert(scip != NULL);
212  	   assert(cutoff != NULL);
213  	   assert(naddedcuts != NULL);
214  	
215  	   if( cutnnz == 0 && SCIPisFeasNegative(scip, cutrhs) ) /*lint !e644*/
216  	   {
217  	      SCIPdebugMsg(scip, " -> gomory cut detected infeasibility with cut 0 <= %g.\n", cutrhs);
218  	      *cutoff = TRUE;
219  	      return SCIP_OKAY;
220  	   }
221  	
222  	   /* Only take efficient cuts, except for cuts with one non-zero coefficient (= bound
223  	    * changes); the latter cuts will be handled internally in sepastore. */
224  	   if( SCIPisEfficacious(scip, cutefficacy) || ( cutnnz == 1 && SCIPisFeasPositive(scip, cutefficacy) ) )
225  	   {
226  	      SCIP_ROW* cut;
227  	      SCIP_SEPA* cutsepa;
228  	      char cutname[SCIP_MAXSTRLEN];
229  	      int v;
230  	
231  	      /* construct cut name */
232  	      if( strongcg )
233  	      {
234  	         cutsepa = sepadata->strongcg;
235  	
236  	         if( c >= 0 )
237  	            (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "scg%" SCIP_LONGINT_FORMAT "_x%d", SCIPgetNLPs(scip), c);
238  	         else
239  	            (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "scg%" SCIP_LONGINT_FORMAT "_s%d", SCIPgetNLPs(scip), -c-1);
240  	      }
241  	      else
242  	      {
243  	         cutsepa = sepadata->gomory;
244  	
245  	         if( c >= 0 )
246  	            (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "gom%" SCIP_LONGINT_FORMAT "_x%d", SCIPgetNLPs(scip), c);
247  	         else
248  	            (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "gom%" SCIP_LONGINT_FORMAT "_s%d", SCIPgetNLPs(scip), -c-1);
249  	      }
250  	
251  	      /* create empty cut */
252  	      SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, cutsepa, cutname, -SCIPinfinity(scip), cutrhs,
253  	            cutislocal, FALSE, sepadata->dynamiccuts) );
254  	
255  	      /* set cut rank */
256  	      SCIProwChgRank(cut, cutrank); /*lint !e644*/
257  	
258  	      /* cache the row extension and only flush them if the cut gets added */
259  	      SCIP_CALL( SCIPcacheRowExtensions(scip, cut) );
260  	
261  	      /* collect all non-zero coefficients */
262  	      for( v = 0; v < cutnnz; ++v )
263  	      {
264  	         SCIP_CALL( SCIPaddVarToRow(scip, cut, vars[cutinds[v]], cutcoefs[v]) );
265  	      }
266  	
267  	      /* flush all changes before adding the cut */
268  	      SCIP_CALL( SCIPflushRowExtensions(scip, cut) );
269  	
270  	      if( SCIProwGetNNonz(cut) == 0 )
271  	      {
272  	         assert( SCIPisFeasNegative(scip, cutrhs) );
273  	         SCIPdebugMsg(scip, " -> gomory cut detected infeasibility with cut 0 <= %g.\n", cutrhs);
274  	         *cutoff = TRUE;
275  	         return SCIP_OKAY;
276  	      }
277  	      else if( SCIProwGetNNonz(cut) == 1 )
278  	      {
279  	         /* Add the bound change as cut to avoid that the LP gets modified. This would mean that the LP is not flushed
280  	          * and the method SCIPgetLPBInvRow() fails; SCIP internally will apply this bound change automatically. */
281  	         SCIP_CALL( SCIPaddRow(scip, cut, TRUE, cutoff) );
282  	         ++(*naddedcuts);
283  	      }
284  	      else
285  	      {
286  	         SCIP_Bool useful;
287  	
288  	         assert(SCIPisInfinity(scip, -SCIProwGetLhs(cut)));
289  	         assert(!SCIPisInfinity(scip, SCIProwGetRhs(cut)));
290  	
291  	         SCIPdebugMsg(scip, " -> %s cut <%s>: rhs=%f, eff=%f\n", strongcg ? "strong-CG" : "gomory", cutname, cutrhs, cutefficacy);
292  	
293  	         SCIP_CALL( evaluateCutNumerics(scip, sepadata, cut, maxdnom, maxscale, &useful) );
294  	
295  	         if( useful )
296  	         {
297  	            SCIPdebugMsg(scip, " -> found %s cut <%s>: act=%f, rhs=%f, norm=%f, eff=%f, min=%f, max=%f (range=%f)\n",
298  	               strongcg ? "strong-CG" : "gomory", cutname, SCIPgetRowLPActivity(scip, cut), SCIProwGetRhs(cut),
299  	               SCIProwGetNorm(cut), SCIPgetCutEfficacy(scip, NULL, cut),
300  	               SCIPgetRowMinCoef(scip, cut), SCIPgetRowMaxCoef(scip, cut),
301  	               SCIPgetRowMaxCoef(scip, cut)/SCIPgetRowMinCoef(scip, cut));
302  	
303  	            if( SCIPisCutNew(scip, cut) )
304  	            {
305  	               /* add global cuts which are not implicit bound changes to the cut pool */
306  	               if( !cutislocal )
307  	               {
308  	                  if( sepadata->delayedcuts )
309  	                  {
310  	                     SCIP_CALL( SCIPaddDelayedPoolCut(scip, cut) );
311  	                  }
312  	                  else
313  	                  {
314  	                     SCIP_CALL( SCIPaddPoolCut(scip, cut) );
315  	                  }
316  	               }
317  	               else
318  	               {
319  	                  /* local cuts we add to the sepastore */
320  	                  SCIP_CALL( SCIPaddRow(scip, cut, FALSE, cutoff) );
321  	               }
322  	
323  	               ++(*naddedcuts);
324  	            }
325  	         }
326  	      }
327  	      /* release the row */
328  	      SCIP_CALL( SCIPreleaseRow(scip, &cut) );
329  	   }
330  	
331  	   return SCIP_OKAY;
332  	}
333  	
334  	/*
335  	 * Callback methods
336  	 */
337  	
338  	/** copy method for separator plugins (called when SCIP copies plugins) */
339  	static
340  	SCIP_DECL_SEPACOPY(sepaCopyGomory)
341  	{  /*lint --e{715}*/
342  	   assert(scip != NULL);
343  	   assert(sepa != NULL);
344  	   assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
345  	
346  	   /* call inclusion method of separator */
347  	   SCIP_CALL( SCIPincludeSepaGomory(scip) );
348  	
349  	   return SCIP_OKAY;
350  	}
351  	
352  	/** destructor of separator to free user data (called when SCIP is exiting) */
353  	/**! [SnippetSepaFreeGomory] */
354  	static
355  	SCIP_DECL_SEPAFREE(sepaFreeGomory)
356  	{  /*lint --e{715}*/
357  	   SCIP_SEPADATA* sepadata;
358  	
359  	   assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
360  	
361  	   /* free separator data */
362  	   sepadata = SCIPsepaGetData(sepa);
363  	   assert(sepadata != NULL);
364  	
365  	   SCIPfreeBlockMemory(scip, &sepadata);
366  	
367  	   SCIPsepaSetData(sepa, NULL);
368  	
369  	   return SCIP_OKAY;
370  	}
371  	/**! [SnippetSepaFreeGomory] */
372  	
373  	/** initialization method of separator (called after problem was transformed) */
374  	static
375  	SCIP_DECL_SEPAINIT(sepaInitGomory)
376  	{
377  	   SCIP_SEPADATA* sepadata;
378  	
379  	   sepadata = SCIPsepaGetData(sepa);
380  	   assert(sepadata != NULL);
381  	
382  	   /* create and initialize random number generator */
383  	   SCIP_CALL( SCIPcreateRandom(scip, &sepadata->randnumgen, DEFAULT_RANDSEED, TRUE) );
384  	
385  	   return SCIP_OKAY;
386  	}
387  	
388  	/** deinitialization method of separator (called before transformed problem is freed) */
389  	static
390  	SCIP_DECL_SEPAEXIT(sepaExitGomory)
391  	{  /*lint --e{715}*/
392  	   SCIP_SEPADATA* sepadata;
393  	
394  	   sepadata = SCIPsepaGetData(sepa);
395  	   assert(sepadata != NULL);
396  	
397  	   SCIPfreeRandom(scip, &sepadata->randnumgen);
398  	
399  	   return SCIP_OKAY;
400  	}
401  	
402  	
403  	/** LP solution separation method of separator */
404  	static
405  	SCIP_DECL_SEPAEXECLP(sepaExeclpGomory)
406  	{  /*lint --e{715}*/
407  	   SCIP_SEPADATA* sepadata;
408  	   SCIP_VAR** vars;
409  	   SCIP_COL** cols;
410  	   SCIP_ROW** rows;
411  	   SCIP_AGGRROW* aggrrow;
412  	   SCIP_VAR* var;
413  	   SCIP_Real* binvrow;
414  	   SCIP_Real* cutcoefs;
415  	   SCIP_Real* basisfrac;
416  	   SCIP_Real* cutefficacies;
417  	   int* basisind;
418  	   int* basisperm;
419  	   int* inds;
420  	   int* cutinds;
421  	   int* colindsproducedcut;
422  	   SCIP_Real maxscale;
423  	   SCIP_Real minfrac;
424  	   SCIP_Real maxfrac;
425  	   SCIP_Real maxcutefficacy;
426  	   SCIP_Longint maxdnom;
427  	   SCIP_Bool cutoff;
428  	   SCIP_Bool separatescg;
429  	   SCIP_Bool separategmi;
430  	   int naddedcuts;
431  	   int nvars;
432  	   int ncols;
433  	   int nrows;
434  	   int ncalls;
435  	   int maxdepth;
436  	   int maxsepacuts;
437  	   int freq;
438  	   int c;
439  	   int i;
440  	   int j;
441  	
442  	   assert(sepa != NULL);
443  	   assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
444  	   assert(scip != NULL);
445  	   assert(result != NULL);
446  	
447  	   *result = SCIP_DIDNOTRUN;
448  	
449  	   sepadata = SCIPsepaGetData(sepa);
450  	   assert(sepadata != NULL);
451  	
452  	   ncalls = SCIPsepaGetNCallsAtNode(sepa);
453  	
454  	   minfrac = sepadata->away;
455  	   maxfrac = 1.0 - sepadata->away;
456  	
457  	   /* only call separator, if we are not close to terminating */
458  	   if( SCIPisStopped(scip) )
459  	      return SCIP_OKAY;
460  	
461  	   /* only call the gomory cut separator a given number of times at each node */
462  	   if( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
463  	      || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
464  	      return SCIP_OKAY;
465  	
466  	   /* only call separator, if an optimal LP solution is at hand */
467  	   if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
468  	      return SCIP_OKAY;
469  	
470  	   /* only call separator, if the LP solution is basic */
471  	   if( !SCIPisLPSolBasic(scip) )
472  	      return SCIP_OKAY;
473  	
474  	   /* only call separator, if there are fractional variables */
475  	   if( SCIPgetNLPBranchCands(scip) == 0 )
476  	      return SCIP_OKAY;
477  	
478  	   /* check whether strong CG cuts should be separated */
479  	   freq = SCIPsepaGetFreq(sepadata->strongcg);
480  	   if( freq > 0 )
481  	      separatescg = (depth % freq == 0);
482  	   else
483  	      separatescg = (freq == depth);
484  	
485  	   /* check whether Gomory MI cuts should be separated */
486  	   freq = SCIPsepaGetFreq(sepadata->gomory);
487  	   if( freq > 0 )
488  	      separategmi = (depth % freq == 0);
489  	   else
490  	      separategmi = (freq == depth);
491  	
492  	   if( !separatescg && !separategmi )
493  	      return SCIP_OKAY;
494  	
495  	   /* get variables data */
496  	   SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
497  	
498  	   /* get LP data */
499  	   SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
500  	   SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
501  	   if( ncols == 0 || nrows == 0 )
502  	      return SCIP_OKAY;
503  	
504  	   /* set the maximal denominator in rational representation of gomory cut and the maximal scale factor to
505  	    * scale resulting cut to integral values to avoid numerical instabilities
506  	    */
507  	   /**@todo find better but still stable gomory cut settings: look at dcmulti, gesa3, khb0525, misc06, p2756 */
508  	   maxdepth = SCIPgetMaxDepth(scip);
509  	   if( depth == 0 )
510  	   {
511  	      maxdnom = 1000;
512  	      maxscale = 1000.0;
513  	   }
514  	   else if( depth <= maxdepth/4 )
515  	   {
516  	      maxdnom = 1000;
517  	      maxscale = 1000.0;
518  	   }
519  	   else if( depth <= maxdepth/2 )
520  	   {
521  	      maxdnom = 100;
522  	      maxscale = 100.0;
523  	   }
524  	   else
525  	   {
526  	      maxdnom = 10;
527  	      maxscale = 10.0;
528  	   }
529  	
530  	   /* allocate temporary memory */
531  	   SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, nvars) );
532  	   SCIP_CALL( SCIPallocBufferArray(scip, &cutinds, nvars) );
533  	   SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) );
534  	   SCIP_CALL( SCIPallocBufferArray(scip, &basisperm, nrows) );
535  	   SCIP_CALL( SCIPallocBufferArray(scip, &basisfrac, nrows) );
536  	   SCIP_CALL( SCIPallocBufferArray(scip, &binvrow, nrows) );
537  	   SCIP_CALL( SCIPallocBufferArray(scip, &inds, nrows) );
538  	   SCIP_CALL( SCIPallocBufferArray(scip, &cutefficacies, nrows) );
539  	   SCIP_CALL( SCIPallocBufferArray(scip, &colindsproducedcut, nrows) );
540  	   SCIP_CALL( SCIPaggrRowCreate(scip, &aggrrow) );
541  	
542  	   /* get basis indices */
543  	   SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) );
544  	
545  	   for( i = 0; i < nrows; ++i )
546  	   {
547  	      SCIP_Real frac = 0.0;
548  	
549  	      c = basisind[i];
550  	      cutefficacies[i] = 0.0;
551  	
552  	      basisperm[i] = i;
553  	
554  	      colindsproducedcut[i] = -1;
555  	
556  	      if( c >= 0 )
557  	      {
558  	
559  	         assert(c < ncols);
560  	         var = SCIPcolGetVar(cols[c]);
561  	         if( SCIPvarGetType(var) != SCIP_VARTYPE_CONTINUOUS )
562  	         {
563  	            frac = SCIPfeasFrac(scip, SCIPcolGetPrimsol(cols[c]));
564  	            frac = MIN(frac, 1.0 - frac);
565  	         }
566  	      }
567  	      else if( sepadata->separaterows )
568  	      {
569  	         SCIP_ROW* row;
570  	
571  	         assert(0 <= -c-1 && -c-1 < nrows);
572  	         row = rows[-c-1];
573  	         if( SCIProwIsIntegral(row) && !SCIProwIsModifiable(row) )
574  	         {
575  	            frac = SCIPfeasFrac(scip, SCIPgetRowActivity(scip, row));
576  	            frac = MIN(frac, 1.0 - frac);
577  	         }
578  	      }
579  	
580  	      if( frac >= minfrac )
581  	      {
582  	         /* slightly change fractionality to have random order for equal fractions */
583  	         basisfrac[i] = frac + SCIPrandomGetReal(sepadata->randnumgen, -1e-6, 1e-6);
584  	      }
585  	      else
586  	      {
587  	         basisfrac[i] = 0.0;
588  	      }
589  	   }
590  	
591  	   /* sort basis indices by fractionality */
592  	   SCIPsortDownRealInt(basisfrac, basisperm, nrows);
593  	
594  	   /* get the maximal number of cuts allowed in a separation round */
595  	   if( depth == 0 )
596  	      maxsepacuts = sepadata->maxsepacutsroot;
597  	   else
598  	      maxsepacuts = sepadata->maxsepacuts;
599  	
600  	   SCIPdebugMsg(scip, "searching gomory cuts: %d cols, %d rows, maxdnom=%" SCIP_LONGINT_FORMAT ", maxscale=%g, maxcuts=%d\n",
601  	      ncols, nrows, maxdnom, maxscale, maxsepacuts);
602  	
603  	   cutoff = FALSE;
604  	   naddedcuts = 0;
605  	
606  	   /* for all basic columns belonging to integer variables, try to generate a gomory cut */
607  	   for( i = 0; i < nrows && naddedcuts < maxsepacuts && !SCIPisStopped(scip) && !cutoff; ++i )
608  	   {
609  	      SCIP_Real cutrhs;
610  	      SCIP_Real cutefficacy = 0.0;
611  	      SCIP_Bool success;
612  	      SCIP_Bool cutislocal;
613  	      SCIP_Bool strongcgsuccess = FALSE;
614  	      int ninds = -1;
615  	      int cutnnz;
616  	      int cutrank;
617  	
618  	      if( basisfrac[i] == 0.0 )
619  	         break;
620  	
621  	      j = basisperm[i];
622  	      c = basisind[j];
623  	
624  	      /* get the row of B^-1 for this basic integer variable with fractional solution value */
625  	      SCIP_CALL( SCIPgetLPBInvRow(scip, j, binvrow, inds, &ninds) );
626  	
627  	      SCIP_CALL( SCIPaggrRowSumRows(scip, aggrrow, binvrow, inds, ninds,
628  	         sepadata->sidetypebasis, allowlocal, 2, (int) MAXAGGRLEN(nvars), &success) );
629  	
630  	      if( !success )
631  	         continue;
632  	
633  	      /* try to create a strong CG cut out of the aggregation row */
634  	      if( separatescg )
635  	      {
636  	         SCIP_CALL( SCIPcalcStrongCG(scip, NULL, POSTPROCESS, BOUNDSWITCH, USEVBDS, allowlocal, minfrac, maxfrac,
637  	            1.0, aggrrow, cutcoefs, &cutrhs, cutinds, &cutnnz, &cutefficacy, &cutrank, &cutislocal, &strongcgsuccess) );
638  	
639  	         /* if we want to generate both cuts, add cut and reset cutefficacy and strongcgsuccess */
640  	         if( strongcgsuccess && sepadata->genbothgomscg )
641  	         {
642  	            assert(allowlocal || !cutislocal); /*lint !e644*/
643  	            SCIP_CALL( addCut(scip, sepadata, vars, c, maxdnom, maxscale, cutnnz, cutinds, cutcoefs, cutefficacy, cutrhs,
644  	                  cutislocal, cutrank, TRUE, &cutoff, &naddedcuts) );
645  	            if( c >= 0 )
646  	            {
647  	               cutefficacies[i] = cutefficacy;
648  	               colindsproducedcut[i] = c;
649  	            }
650  	            cutefficacy = 0.0;
651  	            strongcgsuccess = FALSE;
652  	            if( cutoff )
653  	               break;
654  	         }
655  	      }
656  	
657  	      /* @todo Currently we are using the SCIPcalcMIR() function to compute the coefficients of the Gomory
658  	       *       cut. Alternatively, we could use the direct version (see thesis of Achterberg formula (8.4)) which
659  	       *       leads to cut a of the form \sum a_i x_i \geq 1. Rumor has it that these cuts are better.
660  	       */
661  	
662  	      /* try to create Gomory cut out of the aggregation row */
663  	      if( separategmi )
664  	      {
665  	         /* SCIPcalcMIR will only override the cut if its efficacy is larger than the one of the strongcg cut */
666  	         SCIP_CALL( SCIPcalcMIR(scip, NULL, POSTPROCESS, BOUNDSWITCH, USEVBDS, allowlocal, FIXINTEGRALRHS, NULL, NULL,
667  	            minfrac, maxfrac, 1.0, aggrrow, cutcoefs, &cutrhs, cutinds, &cutnnz, &cutefficacy, &cutrank, &cutislocal, &success) );
668  	
669  	         if( success || strongcgsuccess )
670  	         {
671  	            assert(allowlocal || !cutislocal); /*lint !e644*/
672  	            if( success )
673  	               strongcgsuccess = FALSE;   /* Set strongcgsuccess to FALSE, since the MIR cut has overriden the strongcg cut. */
674  	
675  	            SCIP_CALL( addCut(scip, sepadata, vars, c, maxdnom, maxscale, cutnnz, cutinds, cutcoefs, cutefficacy, cutrhs,
676  	                  cutislocal, cutrank, strongcgsuccess, &cutoff, &naddedcuts) );
677  	            if( c >= 0 )
678  	            {
679  	               cutefficacies[i] = cutefficacy;
680  	               colindsproducedcut[i] = c;
681  	            }
682  	         }
683  	      }
684  	   }
685  	
686  	   /* Add normalized efficacy GMI statistics to history */
687  	   maxcutefficacy = 0.0;
688  	   for( i = 0; i < nrows; ++i )
689  	   {
690  	      if( cutefficacies[i] > maxcutefficacy && colindsproducedcut[i] >= 0 )
691  	      {
692  	         maxcutefficacy = cutefficacies[i];
693  	      }
694  	   }
695  	
696  	
697  	   for( i = 0; i < nrows; ++i )
698  	   {
699  	      if( colindsproducedcut[i] >= 0 && SCIPisEfficacious(scip, cutefficacies[i])  )
700  	      {
701  	         assert( maxcutefficacy > 0.0 );
702  	         var = SCIPcolGetVar(cols[colindsproducedcut[i]]);
703  	         SCIP_CALL( SCIPsetVarLastGMIScore(scip, var, cutefficacies[i] / maxcutefficacy) );
704  	         SCIP_CALL( SCIPincVarGMISumScore(scip, var, cutefficacies[i] / maxcutefficacy) );
705  	      }
706  	   }
707  	
708  	   /* free temporary memory */
709  	   SCIPfreeBufferArray(scip, &inds);
710  	   SCIPfreeBufferArray(scip, &binvrow);
711  	   SCIPfreeBufferArray(scip, &basisfrac);
712  	   SCIPfreeBufferArray(scip, &basisperm);
713  	   SCIPfreeBufferArray(scip, &basisind);
714  	   SCIPfreeBufferArray(scip, &cutinds);
715  	   SCIPfreeBufferArray(scip, &cutcoefs);
716  	   SCIPfreeBufferArray(scip, &cutefficacies);
717  	   SCIPfreeBufferArray(scip, &colindsproducedcut);
718  	   SCIPaggrRowFree(scip, &aggrrow);
719  	
720  	   SCIPdebugMsg(scip, "end searching gomory cuts: found %d cuts\n", naddedcuts);
721  	
722  	   sepadata->lastncutsfound = SCIPgetNCutsFound(scip);
723  	
724  	   /* evaluate the result of the separation */
725  	   if( cutoff )
726  	      *result = SCIP_CUTOFF;
727  	   else if ( naddedcuts > 0 )
728  	      *result = SCIP_SEPARATED;
729  	   else
730  	      *result = SCIP_DIDNOTFIND;
731  	
732  	   return SCIP_OKAY;
733  	}
734  	
735  	
736  	/*
737  	 * separator specific interface methods
738  	 */
739  	
740  	/** LP solution separation method of dummy separator */
741  	static
742  	SCIP_DECL_SEPAEXECLP(sepaExeclpDummy)
743  	{  /*lint --e{715}*/
744  	   assert( result != NULL );
745  	
746  	   *result = SCIP_DIDNOTRUN;
747  	
748  	   return SCIP_OKAY;
749  	}
750  	
751  	/** arbitrary primal solution separation method of dummy separator */
752  	static
753  	SCIP_DECL_SEPAEXECSOL(sepaExecsolDummy)
754  	{  /*lint --e{715}*/
755  	   assert( result != NULL );
756  	
757  	   *result = SCIP_DIDNOTRUN;
758  	
759  	   return SCIP_OKAY;
760  	}
761  	
762  	/** creates the Gomory MIR cut separator and includes it in SCIP */
763  	SCIP_RETCODE SCIPincludeSepaGomory(
764  	   SCIP*                 scip                /**< SCIP data structure */
765  	   )
766  	{
767  	   SCIP_SEPADATA* sepadata;
768  	   SCIP_SEPA* sepa;
769  	
770  	   /* create separator data */
771  	   SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) );
772  	   sepadata->lastncutsfound = 0;
773  	
774  	   /* include separator */
775  	   SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
776  	         SEPA_USESSUBSCIP, SEPA_DELAY,
777  	         sepaExeclpGomory, NULL,
778  	         sepadata) );
779  	   assert(sepa != NULL);
780  	
781  	   SCIP_CALL( SCIPincludeSepaBasic(scip, &sepadata->strongcg, "strongcg", "separator for strong CG cuts", -100000, SEPA_FREQ, 0.0,
782  	      SEPA_USESSUBSCIP, FALSE, sepaExeclpDummy, sepaExecsolDummy, NULL) );
783  	   assert(sepadata->strongcg != NULL);
784  	
785  	   SCIP_CALL( SCIPincludeSepaBasic(scip, &sepadata->gomory, "gomorymi", "separator for Gomory mixed-integer cuts", -100000, SEPA_FREQ, 0.0,
786  	      SEPA_USESSUBSCIP, FALSE, sepaExeclpDummy, sepaExecsolDummy, NULL) );
787  	   assert(sepadata->gomory != NULL);
788  	
789  	   /* set non-NULL pointers to callback methods */
790  	   SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyGomory) );
791  	   SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeGomory) );
792  	   SCIP_CALL( SCIPsetSepaInit(scip, sepa, sepaInitGomory) );
793  	   SCIP_CALL( SCIPsetSepaExit(scip, sepa, sepaExitGomory) );
794  	
795  	   /* mark main separator as a parent */
796  	   SCIPsetSepaIsParentsepa(scip, sepa);
797  	
798  	   /* set pointer from child separators to main separator */
799  	   SCIPsetSepaParentsepa(scip, sepadata->strongcg, sepa);
800  	   SCIPsetSepaParentsepa(scip, sepadata->gomory, sepa);
801  	
802  	   /* add separator parameters */
803  	   SCIP_CALL( SCIPaddIntParam(scip,
804  	         "separating/gomory/maxrounds",
805  	         "maximal number of gomory separation rounds per node (-1: unlimited)",
806  	         &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
807  	   SCIP_CALL( SCIPaddIntParam(scip,
808  	         "separating/gomory/maxroundsroot",
809  	         "maximal number of gomory separation rounds in the root node (-1: unlimited)",
810  	         &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
811  	   SCIP_CALL( SCIPaddIntParam(scip,
812  	         "separating/gomory/maxsepacuts",
813  	         "maximal number of gomory cuts separated per separation round",
814  	         &sepadata->maxsepacuts, FALSE, DEFAULT_MAXSEPACUTS, 0, INT_MAX, NULL, NULL) );
815  	   SCIP_CALL( SCIPaddIntParam(scip,
816  	         "separating/gomory/maxsepacutsroot",
817  	         "maximal number of gomory cuts separated per separation round in the root node",
818  	         &sepadata->maxsepacutsroot, FALSE, DEFAULT_MAXSEPACUTSROOT, 0, INT_MAX, NULL, NULL) );
819  	   SCIP_CALL( SCIPaddIntParam(scip,
820  	         "separating/gomory/maxrank",
821  	         "maximal rank of a gomory cut that could not be scaled to integral coefficients (-1: unlimited)",
822  	         &sepadata->maxrank, FALSE, DEFAULT_MAXRANK, -1, INT_MAX, NULL, NULL) );
823  	   SCIP_CALL( SCIPaddIntParam(scip,
824  	         "separating/gomory/maxrankintegral",
825  	         "maximal rank of a gomory cut that could be scaled to integral coefficients (-1: unlimited)",
826  	         &sepadata->maxrankintegral, FALSE, DEFAULT_MAXRANKINTEGRAL, -1, INT_MAX, NULL, NULL) );
827  	   SCIP_CALL( SCIPaddRealParam(scip,
828  	         "separating/gomory/away",
829  	         "minimal integrality violation of a basis variable in order to try Gomory cut",
830  	         &sepadata->away, FALSE, DEFAULT_AWAY, 1e-4, 0.5, NULL, NULL) );
831  	   SCIP_CALL( SCIPaddBoolParam(scip,
832  	         "separating/gomory/dynamiccuts",
833  	         "should generated cuts be removed from the LP if they are no longer tight?",
834  	         &sepadata->dynamiccuts, FALSE, DEFAULT_DYNAMICCUTS, NULL, NULL) );
835  	   SCIP_CALL( SCIPaddBoolParam(scip,
836  	         "separating/gomory/makeintegral",
837  	         "try to scale cuts to integral coefficients",
838  	         &sepadata->makeintegral, TRUE, DEFAULT_MAKEINTEGRAL, NULL, NULL) );
839  	   SCIP_CALL( SCIPaddBoolParam(scip,
840  	         "separating/gomory/forcecuts",
841  	         "if conversion to integral coefficients failed still consider the cut",
842  	         &sepadata->forcecuts, TRUE, DEFAULT_FORCECUTS, NULL, NULL) );
843  	   SCIP_CALL( SCIPaddBoolParam(scip,
844  	         "separating/gomory/separaterows",
845  	         "separate rows with integral slack",
846  	         &sepadata->separaterows, TRUE, DEFAULT_SEPARATEROWS, NULL, NULL) );
847  	   SCIP_CALL( SCIPaddBoolParam(scip,
848  	         "separating/gomory/delayedcuts",
849  	         "should cuts be added to the delayed cut pool?",
850  	         &sepadata->delayedcuts, TRUE, DEFAULT_DELAYEDCUTS, NULL, NULL) );
851  	   SCIP_CALL( SCIPaddBoolParam(scip,
852  	         "separating/gomory/sidetypebasis",
853  	         "choose side types of row (lhs/rhs) based on basis information?",
854  	         &sepadata->sidetypebasis, TRUE, DEFAULT_SIDETYPEBASIS, NULL, NULL) );
855  	   SCIP_CALL( SCIPaddBoolParam(scip,
856  	         "separating/gomory/trystrongcg",
857  	         "try to generate strengthened Chvatal-Gomory cuts?",
858  	         &sepadata->trystrongcg, TRUE, DEFAULT_TRYSTRONGCG, NULL, NULL) );
859  	   SCIP_CALL( SCIPaddBoolParam(scip,
860  	         "separating/gomory/genbothgomscg",
861  	         "Should both Gomory and strong CG cuts be generated (otherwise take best)?",
862  	         &sepadata->genbothgomscg, TRUE, DEFAULT_GENBOTHGOMSCG, NULL, NULL) );
863  	
864  	   return SCIP_OKAY;
865  	}
866