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_interminor.c
26   	 * @ingroup DEFPLUGINS_SEPA
27   	 * @brief  minor separator with intersection cuts
28   	 * @author Felipe Serrano
29   	 * @author Antonia Chmiela
30   	 *
31   	 * Let X be the matrix of auxiliary variables added for bilinear terms, X_{ij} = x_ix_j.
32   	 * The separator enforces quadratic constraints det(2x2 minor of X) = 0 via intersection cuts.
33   	 */
34   	
35   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
36   	
37   	#include <assert.h>
38   	#include <string.h>
39   	
40   	#include "scip/sepa_interminor.h"
41   	#include "scip/expr.h"
42   	#include "scip/expr_var.h"
43   	#include "scip/expr_pow.h"
44   	#include "scip/expr_product.h"
45   	#include "scip/nlpi_ipopt.h"
46   	#include "scip/cons_nonlinear.h"
47   	
48   	
49   	
50   	#define SEPA_NAME              "interminor"
51   	#define SEPA_DESC              "intersection cuts separator to ensure that 2x2 minors of X (= xx') have determinant 0"
52   	#define SEPA_PRIORITY                 0
53   	#define SEPA_FREQ                    -1
54   	#define SEPA_MAXBOUNDDIST           1.0
55   	#define SEPA_USESSUBSCIP          FALSE /**< does the separator use a secondary SCIP instance? */
56   	#define SEPA_DELAY                FALSE /**< should separation method be delayed, if other separators found cuts? */
57   	
58   	#define DEFAULT_MINCUTVIOL         1e-4 /**< default minimum required violation of a cut */
59   	#define DEFAULT_RANDSEED            157 /**< default random seed */
60   	#define DEFAULT_MAXROUNDS            10 /**< maximal number of separation rounds per node (-1: unlimited) */
61   	#define DEFAULT_MAXROUNDSROOT        -1 /**< maximal number of separation rounds in the root node (-1: unlimited) */
62   	#define BINSEARCH_MAXITERS          120 /**< default iteration limit for binary search */
63   	#define DEFAULT_USESTRENGTHENING  FALSE /**< default for using strengthend intersection cuts to separate */
64   	#define DEFAULT_USEBOUNDS         FALSE /**< default for using nonnegativity bounds when separating */
65   	
66   	/*
67   	 * Data structures
68   	 */
69   	
70   	/** separator data */
71   	struct SCIP_SepaData
72   	{
73   	   SCIP_VAR**            minors;             /**< variables of 2x2 minors; each minor is stored like (auxvar_x^2,auxvar_y^2,auxvar_xy) */
74   	   SCIP_Bool*            isdiagonal;         /**< bool array determining if the variables appearing in the minor are diagonal */
75   	   int                   nminors;            /**< total number of minors */
76   	   int                   minorssize;         /**< size of minors array */
77   	   int                   maxrounds;          /**< maximal number of separation rounds per node (-1: unlimited) */
78   	   int                   maxroundsroot;      /**< maximal number of separation rounds in the root node (-1: unlimited) */
79   	   SCIP_Bool             detectedminors;     /**< has minor detection be called? */
80   	   SCIP_Real             mincutviol;         /**< minimum required violation of a cut */
81   	   SCIP_RANDNUMGEN*      randnumgen;         /**< random number generation */
82   	   SCIP_Bool             usestrengthening;   /**< whether to use strengthened intersection cuts to separate minors */
83   	   SCIP_Bool             usebounds;          /**< whether to also enforce nonegativity bounds of principle minors */
84   	};
85   	
86   	/* these represent a row */
87   	struct rowdata
88   	{
89   	   int*                  vals;               /**< index of the column */
90   	   int                   rowidx;             /**< index corresponding to variable of that row */
91   	   int                   nvals;              /**< number of nonzero entries in column */
92   	   int                   valssize;           /**< size of the array that is currently allocated */
93   	   SCIP_HASHMAP*         auxvars;            /**< entry of the matrix */
94   	};
95   	
96   	/*
97   	 * Local methods
98   	 */
99   	
100  	/** helper method to store a 2x2 minor in the separation data */
101  	static
102  	SCIP_RETCODE sepadataAddMinor(
103  	   SCIP*                 scip,               /**< SCIP data structure */
104  	   SCIP_SEPADATA*        sepadata,           /**< separator data */
105  	   SCIP_VAR*             auxvarxik,          /**< auxiliary variable X_ik = x_i * x_k */
106  	   SCIP_VAR*             auxvarxil,          /**< auxiliary variable X_il = x_i * x_l */
107  	   SCIP_VAR*             auxvarxjk,          /**< auxiliary variable X_jk = x_j * x_k */
108  	   SCIP_VAR*             auxvarxjl,          /**< auxiliary variable X_jl = x_j * x_l */
109  	   SCIP_Bool             isauxvarxikdiag,    /**< is X_ik diagonal? (i.e. i = k) */
110  	   SCIP_Bool             isauxvarxildiag,    /**< is X_il diagonal? (i.e. i = l) */
111  	   SCIP_Bool             isauxvarxjkdiag,    /**< is X_jk diagonal? (i.e. j = k) */
112  	   SCIP_Bool             isauxvarxjldiag     /**< is X_jl diagonal? (i.e. j = l) */
113  	   )
114  	{
115  	   assert(sepadata != NULL);
116  	   assert(auxvarxik != NULL);
117  	   assert(auxvarxil != NULL);
118  	   assert(auxvarxjk != NULL);
119  	   assert(auxvarxjl != NULL);
120  	   assert(auxvarxik != auxvarxil);
121  	   assert(auxvarxjk != auxvarxjl);
122  	
123  	   SCIPdebugMsg(scip, "store 2x2 minor: [%s %s, %s %s]\n", SCIPvarGetName(auxvarxik), SCIPvarGetName(auxvarxil),
124  	         SCIPvarGetName(auxvarxjk), SCIPvarGetName(auxvarxjl));
125  	
126  	   /* reallocate if necessary */
127  	   if( sepadata->minorssize < 4 * (sepadata->nminors + 1) )
128  	   {
129  	      int newsize = SCIPcalcMemGrowSize(scip, 4 * (sepadata->nminors + 1));
130  	      assert(newsize >= 4 * (sepadata->nminors + 1));
131  	
132  	      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(sepadata->minors), sepadata->minorssize, newsize) );
133  	      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(sepadata->isdiagonal), sepadata->minorssize, newsize) );
134  	      sepadata->minorssize = newsize;
135  	   }
136  	
137  	   /* store minor */
138  	   sepadata->minors[4 * sepadata->nminors] = auxvarxik;
139  	   sepadata->minors[4 * sepadata->nminors + 1] = auxvarxil;
140  	   sepadata->minors[4 * sepadata->nminors + 2] = auxvarxjk;
141  	   sepadata->minors[4 * sepadata->nminors + 3] = auxvarxjl;
142  	   sepadata->isdiagonal[4 * sepadata->nminors] = isauxvarxikdiag;
143  	   sepadata->isdiagonal[4 * sepadata->nminors + 1] = isauxvarxildiag;
144  	   sepadata->isdiagonal[4 * sepadata->nminors + 2] = isauxvarxjkdiag;
145  	   sepadata->isdiagonal[4 * sepadata->nminors + 3] = isauxvarxjldiag;
146  	   ++(sepadata->nminors);
147  	
148  	   /* capture variables */
149  	   SCIP_CALL( SCIPcaptureVar(scip, auxvarxik) );
150  	   SCIP_CALL( SCIPcaptureVar(scip, auxvarxil) );
151  	   SCIP_CALL( SCIPcaptureVar(scip, auxvarxjk) );
152  	   SCIP_CALL( SCIPcaptureVar(scip, auxvarxjl) );
153  	
154  	   return SCIP_OKAY;
155  	}
156  	
157  	/** helper method to clear separation data */
158  	static
159  	SCIP_RETCODE sepadataClear(
160  	   SCIP*                 scip,               /**< SCIP data structure */
161  	   SCIP_SEPADATA*        sepadata            /**< separator data */
162  	   )
163  	{
164  	   int i;
165  	
166  	   assert(sepadata != NULL);
167  	
168  	   SCIPdebugMsg(scip, "clear separation data\n");
169  	
170  	   /* release captured variables */
171  	   for( i = 0; i < 4 * sepadata->nminors; ++i )
172  	   {
173  	      assert(sepadata->minors[i] != NULL);
174  	      SCIP_CALL( SCIPreleaseVar(scip, &sepadata->minors[i]) );
175  	   }
176  	
177  	   /* free memory */
178  	   SCIPfreeBlockMemoryArrayNull(scip, &sepadata->minors, sepadata->minorssize);
179  	   SCIPfreeBlockMemoryArrayNull(scip, &sepadata->isdiagonal, sepadata->minorssize);
180  	
181  	   /* reset counters */
182  	   sepadata->nminors = 0;
183  	   sepadata->minorssize = 0;
184  	
185  	   return SCIP_OKAY;
186  	}
187  	
188  	/** helper method to get the variables associated to a minor */
189  	static
190  	SCIP_RETCODE getMinorVars(
191  	   SCIP_SEPADATA*        sepadata,           /**< separator data */
192  	   int                   idx,                /**< index of the stored minor */
193  	   SCIP_VAR**            auxvarxik,          /**< auxiliary variable X_ik = x_i * x_k */
194  	   SCIP_VAR**            auxvarxil,          /**< auxiliary variable X_il = x_i * x_l */
195  	   SCIP_VAR**            auxvarxjk,          /**< auxiliary variable X_jk = x_j * x_k */
196  	   SCIP_VAR**            auxvarxjl,          /**< auxiliary variable X_jl = x_j * x_l */
197  	   SCIP_Bool*            isauxvarxikdiag,    /**< is X_ik diagonal? (i.e. i = k) */
198  	   SCIP_Bool*            isauxvarxildiag,    /**< is X_il diagonal? (i.e. i = l) */
199  	   SCIP_Bool*            isauxvarxjkdiag,    /**< is X_jk diagonal? (i.e. j = k) */
200  	   SCIP_Bool*            isauxvarxjldiag     /**< is X_jl diagonal? (i.e. j = l) */
201  	   )
202  	{
203  	   assert(auxvarxik != NULL);
204  	   assert(auxvarxil != NULL);
205  	   assert(auxvarxjk != NULL);
206  	   assert(auxvarxjl != NULL);
207  	
208  	   *auxvarxik = sepadata->minors[4 * idx];
209  	   *auxvarxil = sepadata->minors[4 * idx + 1];
210  	   *auxvarxjk = sepadata->minors[4 * idx + 2];
211  	   *auxvarxjl = sepadata->minors[4 * idx + 3];
212  	
213  	   *isauxvarxikdiag = sepadata->isdiagonal[4 * idx];
214  	   *isauxvarxildiag = sepadata->isdiagonal[4 * idx + 1];
215  	   *isauxvarxjkdiag = sepadata->isdiagonal[4 * idx + 2];
216  	   *isauxvarxjldiag = sepadata->isdiagonal[4 * idx + 3];
217  	
218  	   return SCIP_OKAY;
219  	}
220  	
221  	
222  	/** adds a new entry (i.e., auxvar) of in (row, col) of matrix M.
223  	 *
224  	 * we have a matrix, M, indexed by the variables
225  	 * M(xi, xk) is the auxiliary variable of xi * xk if it exists
226  	 * We store, for each row of the matrix, the indices of the nonzero column entries (assoc with the given row) and the auxiliary variable for xi * xk
227  	 * The nonzero column entries are stored as an array (struct rowdata)
228  	 * So we have a hashmap mapping each variable (row of the matrix) with its array representing the nonzero entries of the row.
229  	 */
230  	static
231  	SCIP_RETCODE insertIndex(
232  	   SCIP*                 scip,               /**< SCIP data structure */
233  	   SCIP_HASHMAP*         rowmap,             /**< hashmap of the rows of the matrix */
234  	   SCIP_VAR*             row,                /**< variable corresponding to row of new entry */
235  	   SCIP_VAR*             col,                /**< variable corresponding to column of new entry */
236  	   SCIP_VAR*             auxvar,             /**< auxvar to insert into the matrix */
237  	   int*                  rowindices,         /**< array of indices of all variables corresponding to a row */
238  	   int*                  nrows               /**< number of rows */
239  	   )
240  	{
241  	   SCIPdebugMsg(scip, "inserting %s in row %s and col %s \n", SCIPvarGetName(auxvar), SCIPvarGetName(row), SCIPvarGetName(col));
242  	
243  	   /* check whether variable has an array associated to it */
244  	   if( SCIPhashmapExists(rowmap, (void*)row) )
245  	   {
246  	      struct rowdata* arr;
247  	
248  	      arr = (struct rowdata*)SCIPhashmapGetImage(rowmap, (void *)row);
249  	
250  	      /* reallocate if necessary */
251  	      if( arr->valssize < arr->nvals + 1 )
252  	      {
253  	         int newsize = SCIPcalcMemGrowSize(scip, arr->nvals + 1);
254  	         assert(newsize > arr->nvals + 1);
255  	
256  	         SCIP_CALL( SCIPreallocBufferArray(scip, &(arr->vals), newsize) );
257  	         arr->valssize = newsize;
258  	      }
259  	
260  	      /* insert */
261  	      arr->vals[arr->nvals] = SCIPvarGetProbindex(col);
262  	      SCIP_CALL( SCIPhashmapInsert(arr->auxvars, (void*)col, (void *)auxvar) );
263  	      arr->nvals += 1;
264  	   }
265  	   else
266  	   {
267  	      struct rowdata* arr;
268  	
269  	      /* create index array */
270  	      SCIP_CALL( SCIPallocBuffer(scip, &arr) );
271  	      arr->valssize = 10;
272  	      arr->nvals = 0;
273  	      SCIP_CALL( SCIPallocBufferArray(scip, &arr->vals, arr->valssize) );
274  	      SCIP_CALL( SCIPhashmapCreate(&arr->auxvars, SCIPblkmem(scip), arr->valssize) );
275  	
276  	      /* insert */
277  	      arr->rowidx = SCIPvarGetProbindex(row);
278  	      arr->vals[arr->nvals] = SCIPvarGetProbindex(col);
279  	      SCIP_CALL( SCIPhashmapInsert(arr->auxvars, (void*)col, (void *)auxvar) );
280  	      arr->nvals += 1;
281  	
282  	      /* store in hashmap */
283  	      SCIP_CALL( SCIPhashmapInsert(rowmap, (void*)row, (void *)arr) );
284  	
285  	      /* remember the new row */
286  	      rowindices[*nrows] = SCIPvarGetProbindex(row);
287  	      *nrows += 1;
288  	   }
289  	
290  	   return SCIP_OKAY;
291  	}
292  	
293  	/** method to detect and store principal minors */
294  	static
295  	SCIP_RETCODE detectMinors(
296  	   SCIP*                 scip,               /**< SCIP data structure */
297  	   SCIP_SEPADATA*        sepadata            /**< separator data */
298  	   )
299  	{
300  	   SCIP_CONSHDLR* conshdlr;
301  	   SCIP_EXPRITER* it;
302  	   SCIP_HASHMAP* rowmap;
303  	   int* rowvars = NULL;
304  	   int* intersection;
305  	   int nrowvars = 0;
306  	   int c;
307  	   int i;
308  	
309  	#ifdef SCIP_STATISTIC
310  	   SCIP_Real totaltime = -SCIPgetTotalTime(scip);
311  	#endif
312  	
313  	   assert(sepadata != NULL);
314  	
315  	   /* check whether minor detection has been called already */
316  	   if( sepadata->detectedminors )
317  	      return SCIP_OKAY;
318  	
319  	   assert(sepadata->minors == NULL);
320  	   assert(sepadata->nminors == 0);
321  	
322  	   /* we assume that the auxiliary variables in the nonlinear constraint handler have been already generated */
323  	   sepadata->detectedminors = TRUE;
324  	
325  	   /* check whether there are nonlinear constraints available */
326  	   conshdlr = SCIPfindConshdlr(scip, "nonlinear");
327  	   if( conshdlr == NULL || SCIPconshdlrGetNConss(conshdlr) == 0 )
328  	      return SCIP_OKAY;
329  	
330  	   SCIPdebugMsg(scip, "call detectMinors()\n");
331  	
332  	   /* allocate memory */
333  	   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
334  	   SCIP_CALL( SCIPhashmapCreate(&rowmap, SCIPblkmem(scip), SCIPgetNVars(scip)) );
335  	   SCIP_CALL( SCIPallocBufferArray(scip, &rowvars, SCIPgetNVars(scip)) );
336  	   SCIP_CALL( SCIPallocBufferArray(scip, &intersection, SCIPgetNVars(scip)) );
337  	
338  	   /* initialize iterator */
339  	   SCIP_CALL( SCIPexpriterInit(it, NULL, SCIP_EXPRITER_DFS, FALSE) );
340  	
341  	   for( c = 0; c < SCIPconshdlrGetNConss(conshdlr); ++c )
342  	   {
343  	      SCIP_CONS* cons;
344  	      SCIP_EXPR* expr;
345  	      SCIP_EXPR* root;
346  	
347  	      cons = SCIPconshdlrGetConss(conshdlr)[c];
348  	      assert(cons != NULL);
349  	      root = SCIPgetExprNonlinear(cons);
350  	      assert(root != NULL);
351  	
352  	      for( expr = SCIPexpriterRestartDFS(it, root); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) /*lint !e441*//*lint !e440*/
353  	      {
354  	         SCIP_EXPR** children;
355  	         SCIP_VAR* auxvar;
356  	
357  	         SCIPdebugMsg(scip, "visit expression %p in constraint %s\n", (void*)expr, SCIPconsGetName(cons));
358  	
359  	         /* check whether the expression has an auxiliary variable */
360  	         auxvar = SCIPgetExprAuxVarNonlinear(expr);
361  	         if( auxvar == NULL )
362  	         {
363  	            SCIPdebugMsg(scip, "expression has no auxiliary variable -> skip\n");
364  	            continue;
365  	         }
366  	
367  	         children = SCIPexprGetChildren(expr);
368  	
369  	         /* check for expr = (x)^2 */
370  	         if( SCIPexprGetNChildren(expr) == 1 && SCIPisExprPower(scip, expr)
371  	            && SCIPgetExponentExprPow(expr) == 2.0
372  	            && SCIPgetExprAuxVarNonlinear(children[0]) != NULL )
373  	         {
374  	            SCIP_VAR* quadvar;
375  	
376  	            assert(children[0] != NULL);
377  	
378  	            quadvar = SCIPgetExprAuxVarNonlinear(children[0]);
379  	            assert(quadvar != NULL);
380  	
381  	            SCIP_CALL( insertIndex(scip, rowmap, quadvar, quadvar, auxvar, rowvars, &nrowvars) );
382  	         }
383  	         /* check for expr = x_i * x_k */
384  	         else if( SCIPexprGetNChildren(expr) == 2 && SCIPisExprProduct(scip, expr)
385  	            && SCIPgetExprAuxVarNonlinear(children[0]) != NULL && SCIPgetExprAuxVarNonlinear(children[1]) != NULL )
386  	         {
387  	            SCIP_VAR* xi;
388  	            SCIP_VAR* xk;
389  	
390  	            assert(children[0] != NULL);
391  	            assert(children[1] != NULL);
392  	
393  	            xi = SCIPgetExprAuxVarNonlinear(children[0]);
394  	            xk = SCIPgetExprAuxVarNonlinear(children[1]);
395  	
396  	            SCIP_CALL( insertIndex(scip, rowmap, xk, xi, auxvar, rowvars, &nrowvars) );
397  	            SCIP_CALL( insertIndex(scip, rowmap, xi, xk, auxvar, rowvars, &nrowvars) );
398  	         }
399  	      }
400  	   }
401  	
402  	   /* sort the column entries */
403  	   for( i = 0; i < nrowvars; ++i )
404  	   {
405  	      struct rowdata* row;
406  	
407  	      row = (struct rowdata*)SCIPhashmapGetImage(rowmap, (void *)SCIPgetVars(scip)[rowvars[i]]);
408  	      SCIPsortInt(row->vals, row->nvals);
409  	   }
410  	
411  	   /* store 2x2 minors */
412  	   /* TODO: we might store some minors twice since the matrix is symmetric. Handle that! (see unit test for example) */
413  	   for( i = 0; i < nrowvars; ++i )
414  	   {
415  	      int j;
416  	      struct rowdata* rowi;
417  	
418  	      rowi = (struct rowdata*)SCIPhashmapGetImage(rowmap, (void *)SCIPgetVars(scip)[rowvars[i]]);
419  	
420  	      for( j = i + 1; j < nrowvars; ++j )
421  	      {
422  	         struct rowdata* rowj;
423  	         int ninter;
424  	
425  	         rowj = (struct rowdata*)SCIPhashmapGetImage(rowmap, (void *)SCIPgetVars(scip)[rowvars[j]]);
426  	
427  	         SCIPcomputeArraysIntersectionInt(rowi->vals, rowi->nvals, rowj->vals, rowj->nvals, intersection, &ninter);
428  	
429  	         if( ninter > 1)
430  	         {
431  	            int p;
432  	
433  	            for( p = 0; p < ninter - 1; ++p )
434  	            {
435  	               int q;
436  	
437  	               for( q = p + 1; q < ninter; ++q )
438  	               {
439  	                  SCIP_HASHMAP* rowicols;
440  	                  SCIP_HASHMAP* rowjcols;
441  	                  SCIP_VAR* colk;
442  	                  SCIP_VAR* coll;
443  	                  SCIP_VAR* auxvarik;
444  	                  SCIP_VAR* auxvaril;
445  	                  SCIP_VAR* auxvarjk;
446  	                  SCIP_VAR* auxvarjl;
447  	                  int ii;
448  	                  int jj;
449  	                  int k;
450  	                  int l;
451  	                  SCIP_Bool isauxvarikdiag = FALSE;
452  	                  SCIP_Bool isauxvarildiag = FALSE;
453  	                  SCIP_Bool isauxvarjkdiag = FALSE;
454  	                  SCIP_Bool isauxvarjldiag = FALSE;
455  	
456  	                  ii = rowi->rowidx;
457  	                  jj = rowj->rowidx;
458  	                  k = intersection[p];
459  	                  l = intersection[q];
460  	
461  	                  rowicols = rowi->auxvars;
462  	                  rowjcols = rowj->auxvars;
463  	
464  	                  colk = SCIPgetVars(scip)[k];
465  	                  coll = SCIPgetVars(scip)[l];
466  	
467  	                  auxvarik = (SCIP_VAR*) SCIPhashmapGetImage(rowicols, colk);
468  	                  auxvaril = (SCIP_VAR*) SCIPhashmapGetImage(rowicols, coll);
469  	                  auxvarjk = (SCIP_VAR*) SCIPhashmapGetImage(rowjcols, colk);
470  	                  auxvarjl = (SCIP_VAR*) SCIPhashmapGetImage(rowjcols, coll);
471  	
472  	                  if( ii == k )
473  	                     isauxvarikdiag = TRUE;
474  	                  else if( ii == l )
475  	                     isauxvarildiag = TRUE;
476  	                  if( jj == k )
477  	                     isauxvarjkdiag = TRUE;
478  	                  else if( jj == l )
479  	                     isauxvarjldiag = TRUE;
480  	
481  	                  SCIP_CALL( sepadataAddMinor(scip, sepadata, auxvarik, auxvaril, auxvarjk, auxvarjl,
482  	                        isauxvarikdiag, isauxvarildiag, isauxvarjkdiag, isauxvarjldiag) );
483  	               }
484  	            }
485  	         }
486  	      }
487  	      SCIPfreeBufferArrayNull(scip, &rowi->vals);
488  	      SCIPhashmapFree(&rowi->auxvars);
489  	      SCIPfreeBufferArrayNull(scip, &rowi);
490  	   }
491  	
492  	   SCIPdebugMsg(scip, "found %d principal minors in total\n", sepadata->nminors);
493  	
494  	   /* free memory */
495  	   SCIPfreeBufferArray(scip, &intersection);
496  	   SCIPfreeBufferArray(scip, &rowvars);
497  	   SCIPhashmapFree(&rowmap);
498  	   SCIPfreeExpriter(&it);
499  	
500  	#ifdef SCIP_STATISTIC
501  	   totaltime += SCIPgetTotalTime(scip);
502  	   SCIPstatisticMessage("MINOR DETECT %s %f %d %d\n", SCIPgetProbName(scip), totaltime, sepadata->nminors, maxminors);
503  	#endif
504  	
505  	   return SCIP_OKAY;
506  	}
507  	
508  	/** constructs map between lp position of a basic variable and its row in the tableau */
509  	static
510  	SCIP_RETCODE constructBasicVars2TableauRowMap(
511  	   SCIP*                 scip,               /**< SCIP data structure */
512  	   int*                  map                 /**< buffer to store the map */
513  	   )
514  	{
515  	   int* basisind;
516  	   int nrows;
517  	   int i;
518  	
519  	   nrows = SCIPgetNLPRows(scip);
520  	   SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) );
521  	
522  	   SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) );
523  	   for( i = 0; i < nrows; ++i )
524  	   {
525  	      if( basisind[i] >= 0 )
526  	         map[basisind[i]] = i;
527  	   }
528  	
529  	   SCIPfreeBufferArray(scip, &basisind);
530  	
531  	   return SCIP_OKAY;
532  	}
533  	
534  	/** The restriction of the function representing the maximal S-free set to zlp + t * ray has the form
535  	 * SQRT(A t^2 + B t + C) - (D t + E).
536  	 * This function computes the coefficients A, B, C, D, E for the given ray.
537  	 */
538  	static
539  	SCIP_RETCODE computeRestrictionToRay(
540  	   SCIP*                 scip,               /**< SCIP data structure */
541  	   SCIP_Real*            ray,                /**< coefficients of ray */
542  	   SCIP_VAR**            vars,               /**< variables */
543  	   SCIP_Real*            coefs,              /**< buffer to store A, B, C, D, and E of cases 1, 2, 3, or 4a*/
544  	   SCIP_Real*            coefs4b,            /**< buffer to store A, B, C, D, and E of case 4b */
545  	   SCIP_Real*            coefscondition,     /**< buffer to store coefs for checking whether we are in case 4a or 4b */
546  	   SCIP_Bool             usebounds,          /**< TRUE if we want to separate non-negative bound */
547  	   SCIP_Real*            ad,                 /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
548  	   SCIP_Bool*            success             /**< FALSE if we need to abort generation because of numerics */
549  	   )
550  	{
551  	   SCIP_Real eigenvectors[16] = {1.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0};
552  	   SCIP_Real eigenvalues[4] = {0.5, 0.5, -0.5, -0.5};
553  	   SCIP_Real eigencoef = 0.7071067811865475244008443621048490;
554  	   SCIP_Real* a;
555  	   SCIP_Real* b;
556  	   SCIP_Real* c;
557  	   SCIP_Real* d;
558  	   SCIP_Real* e;
559  	   SCIP_Real min;
560  	   SCIP_Real max;
561  	   SCIP_Real norm1;
562  	   SCIP_Real norm2;
563  	   int negidx;
564  	   int posidx;
565  	   int i;
566  	
567  	   *success = TRUE;
568  	
569  	   /* set all coefficients to zero */
570  	   memset(coefs, 0, 5 * sizeof(SCIP_Real));
571  	   memset(coefs4b, 0, 5 * sizeof(SCIP_Real));
572  	   norm1 = 0.0;
573  	   norm2 = 0.0;
574  	
575  	   a = coefs;
576  	   b = coefs + 1;
577  	   c = coefs + 2;
578  	   d = coefs + 3;
579  	   e = coefs + 4;
580  	
581  	   negidx = 2;
582  	   posidx = 0;
583  	   for( i = 0; i < 4; ++i )
584  	   {
585  	      int j;
586  	      SCIP_Real vzlp;
587  	      SCIP_Real vdotray;
588  	
589  	      vzlp = 0;
590  	      vdotray = 0;
591  	
592  	      /* compute eigenvec * ray and eigenvec * solution */
593  	      for( j = 0; j < 4; ++j )
594  	      {
595  	         vdotray += eigencoef * eigenvectors[4 * i + j] * ray[j];
596  	         vzlp += eigencoef * eigenvectors[4 * i + j] * SCIPvarGetLPSol(vars[j]);
597  	      }
598  	
599  	      if( eigenvalues[i] > 0 )
600  	      {
601  	         /* positive eigenvalue: compute D and E */
602  	         *d += eigenvalues[i] * vzlp * vdotray;
603  	         *e += eigenvalues[i] * SQR( vzlp );
604  	
605  	         if( usebounds )
606  	         {
607  	            norm1 += eigenvalues[i] * (1 - SQR( ad[posidx] )) * SQR( vzlp );
608  	            norm2 += SQRT( eigenvalues[i] ) * ad[posidx] * vzlp;
609  	            ++posidx;
610  	         }
611  	      }
612  	      else
613  	      {
614  	         /* negative eigenvalue: compute A, B, and C */
615  	         *a -= eigenvalues[i] * SQR( vdotray );
616  	         *b -= 2.0 * eigenvalues[i] * vzlp * vdotray;
617  	         *c -= eigenvalues[i] * SQR( vzlp );
618  	
619  	         if( usebounds )
620  	         {
621  	            coefs4b[0] -= eigenvalues[i] * (1 - SQR( ad[negidx] )) * SQR( vdotray );
622  	            coefs4b[1] -= 2.0 * eigenvalues[i] * (1 - SQR( ad[negidx] )) * vzlp * vdotray;
623  	            coefs4b[2] -= eigenvalues[i] * (1 - SQR( ad[negidx] )) * SQR( vzlp );
624  	            coefs4b[3] += SQRT( -eigenvalues[i] ) * ad[negidx] * vdotray;
625  	            coefs4b[4] += SQRT( -eigenvalues[i] ) * ad[negidx] * vzlp;
626  	            ++negidx;
627  	         }
628  	      }
629  	   }
630  	
631  	   assert(*e > 0);
632  	
633  	   if( SQRT( *c ) - SQRT( *e ) >= 0.0 )
634  	   {
635  	      assert(SQRT( *c ) - SQRT( *e ) < 1e-6);
636  	      *success = FALSE;
637  	      return SCIP_OKAY;
638  	   }
639  	
640  	   /* finish computation of coefficients when using bounds */
641  	   if( usebounds )
642  	   {
643  	      coefscondition[0] = norm2 / SQRT( *e );
644  	      coefscondition[1] = coefs4b[3];
645  	      coefscondition[2] = coefs4b[4];
646  	
647  	      coefs4b[0] *= norm1 / *e;
648  	      coefs4b[1] *= norm1 / *e;
649  	      coefs4b[2] *= norm1 / *e;
650  	      coefs4b[3] *= norm2 / SQRT( *e );
651  	      coefs4b[4] *= norm2 / SQRT( *e );
652  	
653  	      coefs4b[3] += *d / SQRT( *e );
654  	      coefs4b[4] += SQRT( *e );
655  	
656  	      assert( SQRT( coefs4b[2] ) - coefs4b[4] < 0.0 );
657  	   }
658  	
659  	   /* finish computation of D and E */
660  	   *e = SQRT( *e );
661  	   *d /= *e;
662  	
663  	   /* maybe we want to avoid a large dynamism between A, B and C */
664  	   max = 0.0;
665  	   min = SCIPinfinity(scip);
666  	   for( i = 0; i < 3; ++i )
667  	   {
668  	      SCIP_Real absval;
669  	
670  	      absval = ABS(coefs[i]);
671  	      if( max < absval )
672  	         max = absval;
673  	      if( absval != 0.0 && absval < min )
674  	         min = absval;
675  	   }
676  	
677  	   if( SCIPisHugeValue(scip, max / min) )
678  	   {
679  	#ifdef DEBUG_INTERSECTIONCUT
680  	      printf("Bad numerics: max(A,B,C)/min(A,B,C) is too large (%g)\n", max / min);
681  	#endif
682  	      *success = FALSE;
683  	      return SCIP_OKAY;
684  	   }
685  	
686  	   /* some sanity checks */
687  	   assert(*c >= 0); /* radicand at zero */
688  	   assert(SQRT( *c ) - *e < 0); /* the function at 0 must be negative */
689  	   assert(*a >= 0); /* the function inside the root is convex */
690  	
691  	#ifdef DEBUG_INTERSECTIONCUT
692  	   SCIPinfoMessage(scip, NULL, "Restriction yields: a,b,c,d,e %g %g %g %g %g\n", coefs[0], coefs[1], coefs[2], coefs[3], coefs[4]);
693  	#endif
694  	
695  	   return SCIP_OKAY;
696  	}
697  	
698  	/** returns phi(zlp + t * ray) = SQRT(A t^2 + B t + C) - (D t + E) */  /*lint -e{715}*/
699  	static
700  	SCIP_Real evalPhiAtRay(
701  	   SCIP*                 scip,               /**< SCIP data structure */
702  	   SCIP_Real             t,                  /**< argument of phi restricted to ray */
703  	   SCIP_Real             a,                  /**< value of A */
704  	   SCIP_Real             b,                  /**< value of B */
705  	   SCIP_Real             c,                  /**< value of C */
706  	   SCIP_Real             d,                  /**< value of D */
707  	   SCIP_Real             e                   /**< value of E */
708  	   )
709  	{
710  	#ifdef INTERCUTS_DBLDBL
711  	   SCIP_Real QUAD(lin);
712  	   SCIP_Real QUAD(disc);
713  	   SCIP_Real QUAD(tmp);
714  	   SCIP_Real QUAD(root);
715  	
716  	   /* d * t + e */
717  	   SCIPquadprecProdDD(lin, d, t);
718  	   SCIPquadprecSumQD(lin, lin, e);
719  	
720  	   /* a * t * t */
721  	   SCIPquadprecSquareD(disc, t);
722  	   SCIPquadprecProdQD(disc, disc, a);
723  	
724  	   /* b * t */
725  	   SCIPquadprecProdDD(tmp, b, t);
726  	
727  	   /* a * t * t + b * t */
728  	   SCIPquadprecSumQQ(disc, disc, tmp);
729  	
730  	   /* a * t * t + b * t + c */
731  	   SCIPquadprecSumQD(disc, disc, c);
732  	
733  	   /* sqrt(above): can't take sqrt of 0! */
734  	   if( QUAD_TO_DBL(disc) == 0 )
735  	   {
736  	      QUAD_ASSIGN(root, 0.0);
737  	   }
738  	   else
739  	   {
740  	      SCIPquadprecSqrtQ(root, disc);
741  	   }
742  	
743  	   /* final result */
744  	   QUAD_SCALE(lin, -1.0);
745  	   SCIPquadprecSumQQ(tmp, root, lin);
746  	
747  	   assert(!SCIPisInfinity(scip, t) || QUAD_TO_DBL(tmp) <= 0);
748  	
749  	   return  QUAD_TO_DBL(tmp);
750  	#else
751  	   return SQRT( a * t * t + b * t + c ) - ( d * t + e );
752  	#endif
753  	}
754  	
755  	/** helper function of computeRoot: we want phi to be <= 0 */
756  	static
757  	void doBinarySearch(
758  	   SCIP*                 scip,               /**< SCIP data structure */
759  	   SCIP_Real             a,                  /**< value of A */
760  	   SCIP_Real             b,                  /**< value of B */
761  	   SCIP_Real             c,                  /**< value of C */
762  	   SCIP_Real             d,                  /**< value of D */
763  	   SCIP_Real             e,                  /**< value of E */
764  	   SCIP_Real*            sol                 /**< buffer to store solution; also gives initial point */
765  	   )
766  	{
767  	   SCIP_Real lb = 0.0;
768  	   SCIP_Real ub = *sol;
769  	   SCIP_Real curr;
770  	   int i;
771  	
772  	   for( i = 0; i < BINSEARCH_MAXITERS; ++i )
773  	   {
774  	      SCIP_Real phival;
775  	
776  	      curr = (lb + ub) / 2.0;
777  	      phival = evalPhiAtRay(scip, curr, a, b, c, d, e);
778  	#ifdef INTERCUT_MOREDEBUG
779  	      printf("%d: lb,ub %.10f, %.10f. curr = %g -> phi at curr %g -> phi at lb %g \n", i, lb, ub, curr, phival, evalPhiAtRay(scip, lb, a, b, c, d, e));
780  	#endif
781  	
782  	      if( phival <= 0.0 )
783  	      {
784  	         lb = curr;
785  	         if( SCIPisFeasZero(scip, phival) || SCIPisFeasEQ(scip, ub, lb) )
786  	            break;
787  	      }
788  	      else
789  	         ub = curr;
790  	   }
791  	
792  	   *sol = lb;
793  	}
794  	
795  	/** checks if we are in case 4a, i.e., if
796  	 * (num(xhat_{r+1}(zlp)) / E) * SQRT(A * tsol^2 + B * tsol + C) + w(ray) * tsol + num(yhat_{s+1}(zlp)) <= 0
797  	 */
798  	static
799  	SCIP_Real isCase4a(
800  	   SCIP_Real             tsol,               /**< t in the above formula */
801  	   SCIP_Real*            coefs,              /**< coefficients A, B, C, D, and E of case 4a */
802  	   SCIP_Real*            coefscondition      /**< extra coefficients needed for the evaluation of the condition:
803  	                                              *   num(xhat_{r+1}(zlp)) / E; w(ray); num(yhat_{s+1}(zlp)) */
804  	   )
805  	{
806  	   return (coefscondition[0] * SQRT( coefs[0] * SQR( tsol ) + coefs[1] * tsol + coefs[2] ) + coefscondition[1] *
807  	      tsol + coefscondition[2]) <= 0.0;
808  	}
809  	
810  	/**  finds smallest positive root phi by finding the smallest positive root of
811  	 * (A - D^2) t^2 + (B - 2 D*E) t + (C - E^2) = 0
812  	 *
813  	 * However, we are conservative and want a solution such that phi is negative, but close to 0;
814  	 * thus we correct the result with a binary search
815  	 */
816  	static
817  	SCIP_Real computeRoot(
818  	   SCIP*                 scip,               /**< SCIP data structure */
819  	   SCIP_Real*            coefs               /**< value of A */
820  	   )
821  	{
822  	   SCIP_Real sol;
823  	   SCIP_INTERVAL bounds;
824  	   SCIP_INTERVAL result;
825  	   SCIP_Real a = coefs[0];
826  	   SCIP_Real b = coefs[1];
827  	   SCIP_Real c = coefs[2];
828  	   SCIP_Real d = coefs[3];
829  	   SCIP_Real e = coefs[4];
830  	
831  	   /* there is an intersection point if and only if SQRT(A) > D: here we are beliving in math, this might cause
832  	    * numerical issues
833  	    */
834  	   if( SQRT( a ) <= d )
835  	   {
836  	      sol = SCIPinfinity(scip);
837  	
838  	      return sol;
839  	   }
840  	
841  	   SCIPintervalSetBounds(&bounds, 0.0, SCIPinfinity(scip));
842  	
843  	   /* SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar finds all x such that a x^2 + b x >= c and x in bounds.
844  	    * it is known that if tsol is the root we are looking for, then gamma(zlp + t * ray) <= 0 between 0 and tsol, thus
845  	    * tsol is the smallest t such that (A - D^2) t^2 + (B - 2 D*E) t + (C - E^2) >= 0
846  	    */
847  	   SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_INTERVAL_INFINITY, &result, a - d * d, b - 2.0 * d *
848  	         e, -(c - e * e), bounds);
849  	
850  	   /* it can still be empty because of our infinity, I guess... */
851  	   sol = SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, result) ? SCIPinfinity(scip) : SCIPintervalGetInf(result);
852  	
853  	   /* check that solution is acceptable, ideally it should be <= 0, however when it is positive, we trigger a binary
854  	    * search to make it negative. This binary search might return a solution point that is not at accurately 0 as the
855  	    * one obtained from the function above. Thus, it might fail to satisfy the condition of case 4b in some cases, e.g.,
856  	    * ex8_3_1, bchoco05, etc
857  	    */
858  	   if( evalPhiAtRay(scip, sol, a, b, c, d, e) <= 1e-10 )
859  	   {
860  	#ifdef INTERCUT_MOREDEBUG
861  	      printf("interval solution returned %g -> phival = %g, believe it\n", sol, evalPhiAtRay(sol, a, b, c, d, e));
862  	      printf("don't do bin search\n");
863  	#endif
864  	
865  	      return sol;
866  	   }
867  	   else
868  	   {
869  	      /* perform a binary search to make it negative: this might correct a wrong infinity (e.g. crudeoil_lee1_05) */
870  	#ifdef INTERCUT_MOREDEBUG
871  	      printf("do bin search because phival is %g\n", evalPhiAtRay(scip, sol, a, b, c, d, e));
872  	#endif
873  	      doBinarySearch(scip, a, b, c, d, e, &sol);
874  	   }
875  	
876  	   return sol;
877  	}
878  	
879  	/** The maximal S-free set is gamma(z) <= 0; we find the intersection point of the ray `ray` starting from zlp with the
880  	 * boundary of the S-free set.
881  	 * That is, we find t >= 0 such that gamma(zlp + t * ray) = 0.
882  	 *
883  	 * In cases 1,2, and 3, gamma is of the form
884  	 *    gamma(zlp + t * ray) = SQRT(A t^2 + B t + C) - (D t + E)
885  	 *
886  	 * In the case 4 gamma is of the form
887  	 *    gamma(zlp + t * ray) = SQRT(A t^2 + B t + C) - (D t + E)          if some condition holds
888  	 *                           SQRT(A' t^2 + B' t + C') - (D' t + E')     otherwise
889  	 *
890  	 * It can be shown (given the special properties of gamma) that the smallest positive root of each function of the form
891  	 * SQRT(a t^2 + b t + c) - (d t + e)
892  	 * is the same as the smallest positive root of the quadratic equation:
893  	 *       (SQRT(a t^2 + b t + c) - (d t + e)) * (SQRT(a t^2 + b t + c) + (d t + e)) = 0
894  	 *  <==> (a - d^2) t^2 + (b - 2 d*e) t + (c - e^2) = 0
895  	 *
896  	 * So, in cases 1, 2, and 3, this function just returns the solution of the above equation.
897  	 * In case 4, it first solves the equation assuming we are in the first piece.
898  	 * If there is no solution, then the second piece can't have a solution (first piece >= second piece for all t)
899  	 * Then we check if the solution satisfies the condition.
900  	 * If it doesn't then we solve the equation for the second piece.
901  	 * If it has a solution, then it _has_ to be the solution.
902  	 */
903  	static
904  	SCIP_Real computeIntersectionPoint(
905  	   SCIP*                 scip,               /**< SCIP data structure */
906  	   SCIP_Bool             usebounds,          /**< whether we are in case 4 or not */
907  	   SCIP_Real*            coefs,              /**< values of A, B, C, D, and E of cases 1, 2, 3, or 4a */
908  	   SCIP_Real*            coefs4b,            /**< values of A, B, C, D, and E of case 4b */
909  	   SCIP_Real*            coefscondition      /**< values needed to evaluate condition of case 4 */
910  	   )
911  	{
912  	   SCIP_Real sol;
913  	   SCIP_Real sol4b;
914  	
915  	   assert(coefs != NULL);
916  	
917  	   if( ! usebounds )
918  	      return computeRoot(scip, coefs);
919  	
920  	   assert(coefs4b != NULL);
921  	   assert(coefscondition != NULL);
922  	
923  	   /* compute solution of first piece */
924  	   sol = computeRoot(scip, coefs);
925  	
926  	   /* if there is no solution --> second piece doesn't have solution */
927  	   if( SCIPisInfinity(scip, sol) )
928  	   {
929  	      /* this assert fails on multiplants_mtg5 the problem is that sqrt(A) <= D in 4a but not in 4b,
930  	       * now, this is impossible since the phi4a >= phi4b, so actually sqrt(A) is 10e-15 away from
931  	       * D in 4b
932  	       */
933  	      /* assert(SCIPisInfinity(scip, computeRoot(scip, coefs4b))); */
934  	      return sol;
935  	   }
936  	
937  	   /* if solution of 4a is in 4a, then return */
938  	   if( isCase4a(sol, coefs, coefscondition) )
939  	      return sol;
940  	
941  	   /* not on 4a --> then the intersection point is whatever 4b says: as phi4a >= phi4b, the solution of phi4b should
942  	    * always be larger (but shouldn't be equal at this point given that isCase4a failed, and the condition function
943  	    * evaluates to 0 when phi4a == phi4b) than the solution of phi4a; However, because of numerics (or limits in the
944  	    * binary search) we can find a slightly smaller solution; thus, we just keep the larger one
945  	    */
946  	   sol4b = computeRoot(scip, coefs4b);
947  	
948  	   return MAX(sol, sol4b);
949  	}
950  	
951  	/** adds cutcoef * (col - col*) to rowprep */
952  	static
953  	SCIP_RETCODE addColToCut(
954  	   SCIP*                 scip,               /**< SCIP data structure */
955  	   SCIP_ROWPREP*         rowprep,            /**< rowprep to store intersection cut */
956  	   SCIP_Real             cutcoef,            /**< cut coefficient */
957  	   SCIP_COL*             col                 /**< column to add to rowprep */
958  	   )
959  	{
960  	   assert(col != NULL);
961  	
962  	#ifdef DEBUG_INTERCUTS_NUMERICS
963  	   SCIPinfoMessage(scip, NULL, "adding col %s to cut. %g <= col <= %g\n", SCIPvarGetName(SCIPcolGetVar(col)),
964  	      SCIPvarGetLbLocal(SCIPcolGetVar(col)), SCIPvarGetUbLocal(SCIPcolGetVar(col)));
965  	   SCIPinfoMessage(scip, NULL, "col is active at %s. Value %.15f\n", SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER ? "lower bound" :
966  	      "upper bound" , SCIPcolGetPrimsol(col));
967  	#endif
968  	
969  	   SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPcolGetVar(col), cutcoef) );
970  	   SCIProwprepAddConstant(rowprep, -cutcoef * SCIPcolGetPrimsol(col) );
971  	
972  	   return SCIP_OKAY;
973  	}
974  	
975  	/** adds cutcoef * (slack - slack*) to rowprep
976  	  *
977  	  * row is lhs <= <coefs, vars> + constant <= rhs, thus slack is defined by
978  	  * slack + <coefs, vars> + constant = side
979  	  * If row (slack) is at upper, it means that <coefs,vars*> + constant = rhs, and so
980  	  * slack* = side - rhs --> slack - slack* = rhs - <coefs, vars> - constant.
981  	  * If row (slack) is at lower, then <coefs,vars*> + constant = lhs, and so
982  	  * slack* = side - lhs --> slack - slack* = lhs - <coefs, vars> - constant.
983  	  */
984  	static
985  	SCIP_RETCODE addRowToCut(
986  	   SCIP*                 scip,               /**< SCIP data structure */
987  	   SCIP_ROWPREP*         rowprep,            /**< rowprep to store intersection cut */
988  	   SCIP_Real             cutcoef,            /**< cut coefficient */
989  	   SCIP_ROW*             row,                /**< row, whose slack we are adding to rowprep */
990  	   SCIP_Bool*            success             /**< buffer to store whether the row is nonbasic enough */
991  	   )
992  	{
993  	   int i;
994  	   SCIP_COL** rowcols;
995  	   SCIP_Real* rowcoefs;
996  	   int nnonz;
997  	
998  	   assert(row != NULL);
999  	
1000 	   rowcols = SCIProwGetCols(row);
1001 	   rowcoefs = SCIProwGetVals(row);
1002 	   nnonz = SCIProwGetNLPNonz(row);
1003 	
1004 	#ifdef DEBUG_INTERCUTS_NUMERICS
1005 	   SCIPinfoMessage(scip, NULL, "adding slack var row_%d to cut. %g <= row <= %g\n", SCIProwGetLPPos(row), SCIProwGetLhs(row), SCIProwGetRhs(row));
1006 	   SCIPinfoMessage(scip, NULL, "row is active at %s = %.15f Activity %.15f\n", SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER ? "lhs" :
1007 	   "rhs" , SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER ? SCIProwGetLhs(row) : SCIProwGetRhs(row),
1008 	   SCIPgetRowActivity(scip, row));
1009 	#endif
1010 	
1011 	   if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER )
1012 	   {
1013 	      assert(!SCIPisInfinity(scip, -SCIProwGetLhs(row)));
1014 	      if( ! SCIPisFeasEQ(scip, SCIProwGetLhs(row), SCIPgetRowActivity(scip, row)) )
1015 	      {
1016 	         *success = FALSE;
1017 	         return SCIP_OKAY;
1018 	      }
1019 	
1020 	      SCIProwprepAddConstant(rowprep, SCIProwGetLhs(row) * cutcoef);
1021 	   }
1022 	   else
1023 	   {
1024 	      assert(!SCIPisInfinity(scip, SCIProwGetRhs(row)));
1025 	      if( ! SCIPisFeasEQ(scip, SCIProwGetRhs(row), SCIPgetRowActivity(scip, row)) )
1026 	      {
1027 	         *success = FALSE;
1028 	         return SCIP_OKAY;
1029 	      }
1030 	
1031 	      SCIProwprepAddConstant(rowprep, SCIProwGetRhs(row) * cutcoef);
1032 	   }
1033 	
1034 	   for( i = 0; i < nnonz; i++ )
1035 	   {
1036 	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPcolGetVar(rowcols[i]), -rowcoefs[i] * cutcoef) );
1037 	   }
1038 	
1039 	   SCIProwprepAddConstant(rowprep, -SCIProwGetConstant(row) * cutcoef);
1040 	
1041 	   return SCIP_OKAY;
1042 	}
1043 	
1044 	/** get the tableau rows of the variables in vars */
1045 	static
1046 	SCIP_RETCODE getTableauRows(
1047 	   SCIP*                 scip,               /**< SCIP data structure */
1048 	   SCIP_VAR**            vars,               /**< variables in the minor */
1049 	   int*                  basicvarpos2tableaurow,/**< map between basic var and its tableau row */
1050 	   SCIP_HASHMAP*         tableau,            /**< map between var an its tableau row */
1051 	   SCIP_Real**           tableaurows,        /**< buffer to store tableau row */
1052 	   SCIP_Bool*            success             /**< set to TRUE if no variable had basisstat = ZERO */
1053 	   )
1054 	{
1055 	   int v;
1056 	   int nrows;
1057 	   int ncols;
1058 	
1059 	   *success = TRUE;
1060 	
1061 	   nrows = SCIPgetNLPRows(scip);
1062 	   ncols = SCIPgetNLPCols(scip);
1063 	
1064 	   /* check if we have the tableau row of the variable and if not compute it */
1065 	   for( v = 0; v < 4; ++v )
1066 	   {
1067 	      if( ! SCIPhashmapExists(tableau, (void*)vars[v]) )
1068 	      {
1069 	         SCIP_COL* col;
1070 	
1071 	         /* get column of variable */
1072 	         col = SCIPvarGetCol(vars[v]);
1073 	
1074 	         /* if variable is basic, then get its tableau row and insert it in the hashmap */
1075 	         if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC )
1076 	         {
1077 	            int lppos;
1078 	            SCIP_Real* densetableaurow;
1079 	
1080 	            lppos = SCIPcolGetLPPos(col);
1081 	            SCIP_CALL( SCIPallocBufferArray(scip, &densetableaurow, ncols + nrows) );
1082 	
1083 	            SCIP_CALL( SCIPgetLPBInvRow(scip, basicvarpos2tableaurow[lppos], &densetableaurow[ncols], NULL, NULL) );
1084 	            SCIP_CALL( SCIPgetLPBInvARow(scip, basicvarpos2tableaurow[lppos], &densetableaurow[ncols], densetableaurow, NULL, NULL) );
1085 	
1086 	            /* insert tableau row in hashmap*/
1087 	            SCIP_CALL( SCIPhashmapInsert(tableau, (void*)vars[v], (void *)densetableaurow) );
1088 	         }
1089 	         else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_ZERO )
1090 	         {
1091 	            *success = FALSE;
1092 	            return SCIP_OKAY; /* don't even bother */
1093 	         }
1094 	         else
1095 	         {
1096 	            SCIP_CALL( SCIPhashmapInsert(tableau, (void*)vars[v], (void *)NULL) );
1097 	         }
1098 	      }
1099 	
1100 	      /* get tableau row of var */
1101 	      tableaurows[v] = (SCIP_Real *)SCIPhashmapGetImage(tableau, (void*)vars[v]);
1102 	   }
1103 	   return SCIP_OKAY;
1104 	}
1105 	
1106 	/** computes the cut coefs of the  non-basic (non-slack) variables (correspond to cols) and adds them to the
1107 	 * intersection cut
1108 	 */
1109 	static
1110 	SCIP_RETCODE addCols(
1111 	   SCIP*                 scip,               /**< SCIP data structure */
1112 	   SCIP_VAR**            vars,               /**< variables */
1113 	   SCIP_Real**           tableaurows,        /**< tableau rows corresponding to the variables in vars */
1114 	   SCIP_ROWPREP*         rowprep,            /**< store cut */
1115 	   SCIP_Real*            rays,               /**< buffer to store rays */
1116 	   int*                  nrays,              /**< pointer to store number of nonzero rays */
1117 	   int*                  rayslppos,          /**< buffer to store lppos of nonzero rays */
1118 	   SCIP_Real*            interpoints,        /**< buffer to store intersection points or NULL if not needed */
1119 	   SCIP_Bool             usebounds,          /**< TRUE if we want to separate non-negative bound */
1120 	   SCIP_Real*            ad,                 /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
1121 	   SCIP_Bool*            success             /**< pointer to store whether the generation of cutcoefs was successful */
1122 	   )
1123 	{
1124 	   int i;
1125 	   int ncols;
1126 	   SCIP_COL** cols;
1127 	
1128 	   *success = TRUE;
1129 	
1130 	   /* loop over non-basic (non-slack) variables */
1131 	   cols = SCIPgetLPCols(scip);
1132 	   ncols = SCIPgetNLPCols(scip);
1133 	   for( i = 0; i < ncols; ++i )
1134 	   {
1135 	      SCIP_COL* col;
1136 	      SCIP_Real coefs[5];
1137 	      SCIP_Real coefs4b[5];
1138 	      SCIP_Real coefscondition[3];
1139 	      SCIP_Real factor;
1140 	      SCIP_Bool israynonzero;
1141 	      SCIP_Real cutcoef;
1142 	      SCIP_Real interpoint;
1143 	      int v;
1144 	
1145 	      col = cols[i];
1146 	
1147 	      /* set factor to store entries of ray as = [-BinvL, BinvU] */
1148 	      if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER )
1149 	         factor = -1.0;
1150 	      else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER )
1151 	         factor = 1.0;
1152 	      else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_ZERO )
1153 	      {
1154 	         *success = FALSE;
1155 	         return SCIP_OKAY;
1156 	      }
1157 	      else
1158 	         continue;
1159 	
1160 	      /* build the ray */
1161 	      israynonzero = FALSE;
1162 	      for( v = 0; v < 4; ++v )
1163 	      {
1164 	         if( tableaurows[v] != NULL )
1165 	            rays[(*nrays) * 4 + v] = factor * (SCIPisZero(scip, tableaurows[v][i]) ? 0.0 : tableaurows[v][i]);
1166 	         else
1167 	         {
1168 	            if( col == SCIPvarGetCol(vars[v]) )
1169 	               rays[(*nrays) * 4 + v] = -factor;
1170 	            else
1171 	               rays[(*nrays) * 4 + v] = 0.0;
1172 	         }
1173 	
1174 	         israynonzero = israynonzero || (rays[(*nrays) * 4 + v] != 0.0);
1175 	      }
1176 	
1177 	      /* do nothing if ray is 0 */
1178 	      if( ! israynonzero )
1179 	         continue;
1180 	
1181 	      /* compute the cut */
1182 	      SCIP_CALL( computeRestrictionToRay(scip, &rays[(*nrays) * 4], vars, coefs, coefs4b, coefscondition, usebounds,
1183 	            ad, success) );
1184 	
1185 	      if( *success == FALSE )
1186 	         return SCIP_OKAY;
1187 	
1188 	      /* compute intersection point */
1189 	      interpoint = computeIntersectionPoint(scip, usebounds, coefs, coefs4b, coefscondition);
1190 	
1191 	      /* store intersection points */
1192 	      interpoints[*nrays] = interpoint;
1193 	
1194 	      /* remember lppos */
1195 	      rayslppos[*nrays] = i;
1196 	
1197 	      /* count nonzero rays */
1198 	      *nrays += 1;
1199 	
1200 	      /* compute cut coef */
1201 	      cutcoef = SCIPisInfinity(scip, interpoint) ? 0.0 : 1.0 / interpoint;
1202 	
1203 	      /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
1204 	      assert(SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER || SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER);
1205 	      SCIP_CALL( addColToCut(scip, rowprep, SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER ? -cutcoef :
1206 	            cutcoef, col) );
1207 	   }
1208 	
1209 	   return SCIP_OKAY;
1210 	}
1211 	
1212 	/** computes the cut coefs of the non-basic slack variables (correspond to rows) and adds them to the
1213 	 * intersection cut
1214 	 */
1215 	static
1216 	SCIP_RETCODE addRows(
1217 	   SCIP*                 scip,               /**< SCIP data structure */
1218 	   SCIP_VAR**            vars,               /**< variables */
1219 	   SCIP_Real**           tableaurows,        /**< tableau rows corresponding to the variables in vars */
1220 	   SCIP_ROWPREP*         rowprep,            /**< store cut */
1221 	   SCIP_Real*            rays,               /**< buffer to store rays */
1222 	   int*                  nrays,              /**< pointer to store number of nonzero rays */
1223 	   int*                  rayslppos,          /**< buffer to store lppos of nonzero rays */
1224 	   SCIP_Real*            interpoints,        /**< buffer to store intersection points or NULL if not needed */
1225 	   SCIP_Bool             usebounds,          /**< TRUE if we want to separate non-negative bound */
1226 	   SCIP_Real*            ad,                 /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
1227 	   SCIP_Bool*            success             /**< pointer to store whether the generation of cutcoefs was successful */
1228 	   )
1229 	{
1230 	   int i;
1231 	   int nrows;
1232 	   int ncols;
1233 	   SCIP_ROW** rows;
1234 	
1235 	   nrows = SCIPgetNLPRows(scip);
1236 	   ncols = SCIPgetNLPCols(scip);
1237 	
1238 	   *success = TRUE;
1239 	
1240 	   /* loop over non-basic slack variables */
1241 	   rows = SCIPgetLPRows(scip);
1242 	   for( i = 0; i < nrows; ++i )
1243 	   {
1244 	      SCIP_ROW* row;
1245 	      SCIP_Real coefs[5];
1246 	      SCIP_Real coefs4b[5];
1247 	      SCIP_Real coefscondition[3];
1248 	      SCIP_Real factor;
1249 	      SCIP_Bool israynonzero;
1250 	      SCIP_Real cutcoef;
1251 	      SCIP_Real interpoint;
1252 	      int v;
1253 	
1254 	      row = rows[i];
1255 	
1256 	      /* set factor to store entries of ray as = [BinvL, -BinvU] */
1257 	      if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER )
1258 	         factor = 1.0;
1259 	      else if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER )
1260 	         factor = -1.0;
1261 	      else if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_ZERO )
1262 	      {
1263 	         *success = FALSE;
1264 	         return SCIP_OKAY;
1265 	      }
1266 	      else
1267 	         continue;
1268 	
1269 	      /* build the ray */
1270 	      israynonzero = FALSE;
1271 	      for( v = 0; v < 4; ++v )
1272 	      {
1273 	         int idx;
1274 	
1275 	         idx = ncols + i;
1276 	
1277 	         if( tableaurows[v] != NULL )
1278 	            rays[(*nrays) * 4 + v] = factor * (SCIPisZero(scip, tableaurows[v][idx]) ? 0.0 : tableaurows[v][idx]);
1279 	         else
1280 	         {
1281 	            /* TODO: We assume that slack variables can never occure in the minor. This is correct, right? */
1282 	            rays[(*nrays) * 4 + v] = 0.0;
1283 	         }
1284 	
1285 	         israynonzero = israynonzero || (rays[(*nrays) * 4 + v] != 0.0);
1286 	      }
1287 	
1288 	      /* do nothing if ray is 0 */
1289 	      if( ! israynonzero )
1290 	         continue;
1291 	
1292 	      /* compute the cut */
1293 	      SCIP_CALL( computeRestrictionToRay(scip, &rays[(*nrays) * 4], vars, coefs, coefs4b, coefscondition, usebounds,
1294 	            ad, success) );
1295 	
1296 	      if( *success == FALSE )
1297 	         return SCIP_OKAY;
1298 	
1299 	      /* compute intersection point */
1300 	      interpoint = computeIntersectionPoint(scip, usebounds, coefs, coefs4b, coefscondition);
1301 	
1302 	      /* store intersection points */
1303 	      interpoints[*nrays] = interpoint;
1304 	
1305 	      /* store lppos of ray, make it negative so we can differentiate between cols and rows */
1306 	      rayslppos[*nrays] = -i - 1;
1307 	
1308 	      /* count nonzero rays */
1309 	      *nrays += 1;
1310 	
1311 	      /* compute cut coef */
1312 	      cutcoef = SCIPisInfinity(scip, interpoint) ? 0.0 : 1.0 / interpoint;
1313 	
1314 	      /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
1315 	      assert(SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER);
1316 	
1317 	      SCIP_CALL( addRowToCut(scip, rowprep, SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER ? cutcoef :
1318 	            -cutcoef, row, success) ); /* rows have flipper base status! */
1319 	   }
1320 	
1321 	   return SCIP_OKAY;
1322 	}
1323 	
1324 	/* checks if two rays are linearly dependent */
1325 	static
1326 	SCIP_Bool raysAreDependent(
1327 	   SCIP*                 scip,               /**< SCIP data structure */
1328 	   SCIP_Real*            ray1,               /**< coefficients of ray 1 */
1329 	   SCIP_Real*            ray2,               /**< coefficients of ray 2 */
1330 	   SCIP_Real*            coef                /**< pointer to store coef (s.t. r1 = coef * r2) in case rays are dependent */
1331 	   )
1332 	{
1333 	   int i;
1334 	
1335 	   *coef = 0.0;
1336 	
1337 	   for( i = 0; i < 4; ++i )
1338 	   {
1339 	      /* rays cannot be dependent if one ray has zero entry and the other one doesn't */
1340 	      if( (SCIPisZero(scip, ray1[i]) && ! SCIPisZero(scip, ray2[i])) ||
1341 	         (! SCIPisZero(scip, ray1[i]) && SCIPisZero(scip, ray2[i])) )
1342 	      {
1343 	         return FALSE;
1344 	      }
1345 	
1346 	      if( *coef != 0.0 )
1347 	      {
1348 	         /* cannot be dependent if the coefs aren't equal for all entries */
1349 	         if( ! SCIPisFeasEQ(scip, *coef, ray1[i] / ray2[i]) )
1350 	            return FALSE;
1351 	      }
1352 	      else
1353 	         *coef = ray1[i] / ray2[i];
1354 	   }
1355 	
1356 	   return TRUE;
1357 	}
1358 	
1359 	/** finds the smallest negative steplength for the current ray r_idx such that the combination
1360 	 * of r_idx with all rays not in the recession cone is in the recession cone
1361 	 */
1362 	static
1363 	SCIP_RETCODE findRho(
1364 	   SCIP*                 scip,               /**< SCIP data structure */
1365 	   SCIP_Real*            rays,               /**< rays */
1366 	   int                   nrays,              /**< number of nonzero rays */
1367 	   int                   idx,                /**< index of current ray we want to find rho for */
1368 	   SCIP_Real*            interpoints,        /**< intersection points of nonzero rays */
1369 	   SCIP_VAR**            vars,               /**< variables */
1370 	   SCIP_Real*            rho,                /**< pointer to store the optimal rho */
1371 	   SCIP_Bool             usebounds,          /**< TRUE if we want to separate non-negative bound */
1372 	   SCIP_Real*            ad,                 /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
1373 	   SCIP_Bool*            success             /**< TRUE if computation of rho was successful */
1374 	   )
1375 	{
1376 	   int i;
1377 	
1378 	   *success = TRUE;
1379 	
1380 	   /* go through all rays not in the recession cone and compute the largest negative steplength possible. The
1381 	    * smallest of them is then the steplength rho we use for the current ray */
1382 	   *rho = 0;
1383 	   for( i = 0; i < nrays; ++i )
1384 	   {
1385 	      SCIP_Real currentrho;
1386 	      SCIP_Real coef;
1387 	
1388 	      if( SCIPisInfinity(scip, interpoints[i]) )
1389 	         continue;
1390 	
1391 	      /* if the rays are linearly independent, we don't need to search for rho */
1392 	      if( raysAreDependent(scip, &rays[4 * i], &rays[4 * idx], &coef) )
1393 	         currentrho = coef * interpoints[i];
1394 	      else
1395 	      {
1396 	         SCIP_Real lb;
1397 	         SCIP_Real ub;
1398 	         SCIP_Real alpha;
1399 	         int j;
1400 	
1401 	         /* do binary search by lookig at the convex combinations of r_i and r_j */
1402 	         lb = 0.0;
1403 	         ub = 1.0;
1404 	
1405 	         for( j = 0; j < BINSEARCH_MAXITERS; ++j )
1406 	         {
1407 	            SCIP_Real coefs[5];
1408 	            SCIP_Real coefs4b[5];
1409 	            SCIP_Real coefscondition[3];
1410 	            SCIP_Real newray[4];
1411 	            SCIP_Real interpoint;
1412 	            int k;
1413 	
1414 	            alpha = (lb + ub) / 2.0;
1415 	
1416 	            /* build the ray alpha * ray_i + (1 - alpha) * ray_idx */
1417 	            for( k = 0; k < 4; ++k )
1418 	               newray[k] = alpha * rays[4 * i + k] - (1 - alpha) * rays[4 * idx + k];
1419 	
1420 	            /* restrict phi to the "new" ray */
1421 	            SCIP_CALL( computeRestrictionToRay(scip, newray, vars, coefs, coefs4b, coefscondition, usebounds,
1422 	                  ad, success) );
1423 	
1424 	            if( ! *success )
1425 	               return SCIP_OKAY;
1426 	
1427 	            /* check if restriction to "new" ray is numerically nasty. If so, treat the corresponding rho as if phi is
1428 	             * positive
1429 	             */
1430 	
1431 	            /* compute intersection point */
1432 	            interpoint = computeIntersectionPoint(scip, usebounds, coefs, coefs4b, coefscondition);
1433 	
1434 	            /* no root exists */
1435 	            if( SCIPisInfinity(scip, interpoint) )
1436 	            {
1437 	               lb = alpha;
1438 	               if( SCIPisEQ(scip, ub, lb) )
1439 	                  break;
1440 	            }
1441 	            else
1442 	               ub = alpha;
1443 	         }
1444 	
1445 	         /* now we found the best convex combination which we use to derive the corresponding coef. If alpha = 0, we
1446 	          * cannot move the ray in the recession cone, i.e. strengthening is not possible */
1447 	         if( SCIPisZero(scip, alpha) )
1448 	         {
1449 	            *rho = -SCIPinfinity(scip);
1450 	            return SCIP_OKAY;
1451 	         }
1452 	         else
1453 	            currentrho = (alpha - 1) * interpoints[i] / alpha;
1454 	      }
1455 	
1456 	      if( currentrho < *rho )
1457 	         *rho = currentrho;
1458 	   }
1459 	
1460 	   return SCIP_OKAY;
1461 	}
1462 	
1463 	/** computes negative steplengths for the rays that are in the recession cone of the S-free set, i.e.,
1464 	 * which have an infinite intersection point.
1465 	 */
1466 	static
1467 	SCIP_RETCODE computeNegCutcoefs(
1468 	   SCIP*                 scip,               /**< SCIP data structure */
1469 	   SCIP_VAR**            vars,               /**< variables */
1470 	   SCIP_Real*            rays,               /**< rays */
1471 	   int                   nrays,              /**< number of nonzero rays */
1472 	   int*                  rayslppos,          /**< lppos of nonzero rays */
1473 	   SCIP_Real*            interpoints,        /**< intersection points */
1474 	   SCIP_ROWPREP*         rowprep,            /**< rowprep for the generated cut */
1475 	   SCIP_Bool             usebounds,          /**< TRUE if we want to separate non-negative bound */
1476 	   SCIP_Real*            ad,                 /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
1477 	   SCIP_Bool*            success             /**< if a cut candidate could be computed */
1478 	   )
1479 	{
1480 	   SCIP_COL** cols;
1481 	   SCIP_ROW** rows;
1482 	   int i;
1483 	
1484 	   *success = TRUE;
1485 	
1486 	   cols = SCIPgetLPCols(scip);
1487 	   rows = SCIPgetLPRows(scip);
1488 	
1489 	   /* go through all intersection points that are equal to infinity -> these correspond to the rays which are in the
1490 	    * recession cone of C, i.e. the rays for which we (possibly) can compute a negative steplength */
1491 	   for( i = 0; i < nrays ; ++i )
1492 	   {
1493 	      SCIP_Real rho;
1494 	      SCIP_Real cutcoef;
1495 	      int lppos;
1496 	
1497 	      if( !SCIPisInfinity(scip, interpoints[i]) )
1498 	         continue;
1499 	
1500 	      /* compute the smallest rho */
1501 	      SCIP_CALL( findRho(scip, rays, nrays, i, interpoints, vars, &rho, usebounds, ad, success) );
1502 	
1503 	      if( ! *success )
1504 	         continue;
1505 	
1506 	      /* compute cut coef */
1507 	      cutcoef = SCIPisInfinity(scip, -rho) ? 0.0 : 1.0 / rho;
1508 	
1509 	      /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
1510 	      lppos = rayslppos[i];
1511 	      if( lppos < 0 )
1512 	      {
1513 	         lppos = -lppos - 1;
1514 	
1515 	         assert(SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(rows[lppos]) ==
1516 	               SCIP_BASESTAT_UPPER);
1517 	
1518 	         SCIP_CALL( addRowToCut(scip, rowprep, SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_UPPER ? cutcoef :
1519 	                  -cutcoef, rows[lppos], success) ); /* rows have flipped base status! */
1520 	
1521 	         if( ! *success )
1522 	            return SCIP_OKAY;
1523 	      }
1524 	      else
1525 	      {
1526 	         assert(SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER || SCIPcolGetBasisStatus(cols[lppos]) ==
1527 	               SCIP_BASESTAT_LOWER);
1528 	         SCIP_CALL( addColToCut(scip, rowprep, SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER ? -cutcoef :
1529 	                  cutcoef, cols[lppos]) );
1530 	      }
1531 	   }
1532 	
1533 	   return SCIP_OKAY;
1534 	}
1535 	
1536 	/** separates cuts for stored principal minors */
1537 	static
1538 	SCIP_RETCODE separateDeterminant(
1539 	   SCIP*                 scip,               /**< SCIP data structure */
1540 	   SCIP_SEPA*            sepa,               /**< separator */
1541 	   SCIP_SEPADATA*        sepadata,           /**< separator data */
1542 	   SCIP_VAR*             xik,                /**< variable X_ik = x_i * x_k */
1543 	   SCIP_VAR*             xil,                /**< variable X_il = x_i * x_l */
1544 	   SCIP_VAR*             xjk,                /**< variable X_jk = x_j * x_k */
1545 	   SCIP_VAR*             xjl,                /**< variable X_jl = x_j * x_l */
1546 	   SCIP_Bool*            isxikdiag,          /**< is X_ik diagonal? (i.e. i = k) */
1547 	   SCIP_Bool*            isxildiag,          /**< is X_il diagonal? (i.e. i = l) */
1548 	   SCIP_Bool*            isxjkdiag,          /**< is X_jk diagonal? (i.e. j = k) */
1549 	   SCIP_Bool*            isxjldiag,          /**< is X_jl diagonal? (i.e. j = l) */
1550 	   int*                  basicvarpos2tableaurow,/**< map from basic var to its tableau row */
1551 	   SCIP_HASHMAP*         tableau,            /**< map from var to its tableau row */
1552 	   SCIP_RESULT*          result              /**< pointer to store the result of the separation call */
1553 	   )
1554 	{
1555 	   SCIP_ROWPREP* rowprep;
1556 	   SCIP_VAR* vars[4] = {xik, xjl, xil, xjk};
1557 	   SCIP_Real* tableaurows[4];
1558 	   SCIP_Real* interpoints;
1559 	   SCIP_Real* rays;
1560 	   int nrays;
1561 	   int* rayslppos;
1562 	   int ncols;
1563 	   int nrows;
1564 	   SCIP_Bool success;
1565 	   SCIP_Real ad[4] = {0.0, 0.0, 0.0, 0.0};
1566 	   SCIP_Real solxik;
1567 	   SCIP_Real solxil;
1568 	   SCIP_Real solxjk;
1569 	   SCIP_Real solxjl;
1570 	
1571 	   ncols = SCIPgetNLPCols(scip);
1572 	   nrows = SCIPgetNLPRows(scip);
1573 	
1574 	   /* allocate memory for intersection points */
1575 	   SCIP_CALL( SCIPallocBufferArray(scip, &interpoints, ncols + nrows) );
1576 	
1577 	   /* allocate memory for rays */
1578 	   SCIP_CALL( SCIPallocBufferArray(scip, &rays, 4 * (ncols + nrows)) );
1579 	   SCIP_CALL( SCIPallocBufferArray(scip, &rayslppos, ncols + nrows) );
1580 	
1581 	   /* cut (in the nonbasic space) is of the form alpha^T x >= 1 */
1582 	   SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, SCIP_SIDETYPE_LEFT, TRUE) );
1583 	   SCIProwprepAddSide(rowprep, 1.0);
1584 	
1585 	   /* check if we have the tableau row of the variable and if not compute it */
1586 	   SCIP_CALL( getTableauRows(scip, vars, basicvarpos2tableaurow, tableau, tableaurows, &success) );
1587 	
1588 	   if( ! success )
1589 	      goto CLEANUP;
1590 	
1591 	   /* if we want to enforce bounds, set the right a and d to enforce aTx + dTy <= 0 */
1592 	   if( sepadata->usebounds )
1593 	   {
1594 	      solxik = SCIPvarGetLPSol(xik);
1595 	      solxil = SCIPvarGetLPSol(xil);
1596 	      solxjk = SCIPvarGetLPSol(xjk);
1597 	      solxjl = SCIPvarGetLPSol(xjl);
1598 	
1599 	      if( isxikdiag && SCIPisFeasNegative(scip, solxik) )
1600 	      {
1601 	         ad[0] = -1.0;
1602 	         ad[2] = 1.0;
1603 	      }
1604 	      else if( isxjldiag && SCIPisFeasNegative(scip, solxjl) )
1605 	      {
1606 	         ad[0] = -1.0;
1607 	         ad[2] = -1.0;
1608 	      }
1609 	      else if( isxildiag && SCIPisFeasNegative(scip, solxil) )
1610 	      {
1611 	         ad[1] = 1.0;
1612 	         ad[3] = -1.0;
1613 	      }
1614 	      else if( isxjkdiag && SCIPisFeasNegative(scip, solxjk) )
1615 	      {
1616 	         ad[1] = -1.0;
1617 	         ad[3] = -1.0;
1618 	      }
1619 	   }
1620 	
1621 	   nrays = 0;
1622 	   /* loop over each non-basic var; get the ray; compute cut coefficient */
1623 	   SCIP_CALL( addCols(scip, vars, tableaurows, rowprep, rays, &nrays, rayslppos, interpoints, sepadata->usebounds, ad, &success) );
1624 	
1625 	   if( ! success )
1626 	      goto CLEANUP;
1627 	
1628 	   /* loop over non-basic slack variables */
1629 	   SCIP_CALL( addRows(scip, vars, tableaurows, rowprep, rays, &nrays, rayslppos, interpoints, sepadata->usebounds, ad, &success) );
1630 	
1631 	   if( ! success )
1632 	      goto CLEANUP;
1633 	
1634 	   /* do strengthening */
1635 	   if( sepadata->usestrengthening )
1636 	   {
1637 	      SCIP_CALL( computeNegCutcoefs(scip, vars, rays, nrays, rayslppos, interpoints, rowprep, sepadata->usebounds, ad, &success) );
1638 	
1639 	      if( ! success )
1640 	         goto CLEANUP;
1641 	   }
1642 	
1643 	   /* merge coefficients that belong to same variable */
1644 	   SCIPmergeRowprepTerms(scip, rowprep);
1645 	
1646 	   SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, NULL, sepadata->mincutviol, NULL, &success) );
1647 	
1648 	   /* if cleanup was successfull, create row out of rowprep and add it */
1649 	   if( success )
1650 	   {
1651 	      SCIP_ROW* row;
1652 	      SCIP_Bool infeasible;
1653 	
1654 	      /* create row */
1655 	      SCIP_CALL( SCIPgetRowprepRowSepa(scip, &row, rowprep, sepa) );
1656 	
1657 	      assert(SCIPgetCutEfficacy(scip, NULL, row) > 0.0);
1658 	
1659 	      /* add row */
1660 	      SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
1661 	
1662 	      if( infeasible )
1663 	         *result = SCIP_CUTOFF;
1664 	      else
1665 	         *result = SCIP_SEPARATED;
1666 	
1667 	      SCIP_CALL( SCIPreleaseRow(scip, &row) );
1668 	   }
1669 	
1670 	CLEANUP:
1671 	   SCIPfreeRowprep(scip, &rowprep);
1672 	   SCIPfreeBuffer(scip, &rayslppos);
1673 	   SCIPfreeBuffer(scip, &rays);
1674 	   SCIPfreeBuffer(scip, &interpoints);
1675 	
1676 	   return SCIP_OKAY;
1677 	}
1678 	
1679 	
1680 	/** separates cuts for stored principal minors */
1681 	static
1682 	SCIP_RETCODE separatePoint(
1683 	   SCIP*                 scip,               /**< SCIP data structure */
1684 	   SCIP_SEPA*            sepa,               /**< separator */
1685 	   SCIP_RESULT*          result              /**< pointer to store the result of the separation call */
1686 	   )
1687 	{
1688 	   SCIP_SEPADATA* sepadata;
1689 	   SCIP_HASHMAP* tableau = NULL;
1690 	   int* basicvarpos2tableaurow = NULL; /* map between basic var and its tableau row */
1691 	   int i;
1692 	
1693 	   assert(sepa != NULL);
1694 	   assert(result != NULL);
1695 	
1696 	   *result = SCIP_DIDNOTRUN;
1697 	
1698 	   sepadata = SCIPsepaGetData(sepa);
1699 	   assert(sepadata != NULL);
1700 	
1701 	   /* check whether there are some minors available */
1702 	   if( sepadata->nminors == 0 )
1703 	      return SCIP_OKAY;
1704 	
1705 	   *result = SCIP_DIDNOTFIND;
1706 	
1707 	   /* loop over the minors and if they are violated build cut */
1708 	   for( i = 0; i < sepadata->nminors && (*result != SCIP_CUTOFF); ++i )
1709 	   {
1710 	      SCIP_VAR* auxvarxik;
1711 	      SCIP_VAR* auxvarxil;
1712 	      SCIP_VAR* auxvarxjk;
1713 	      SCIP_VAR* auxvarxjl;
1714 	      SCIP_Bool isauxvarxikdiag;
1715 	      SCIP_Bool isauxvarxildiag;
1716 	      SCIP_Bool isauxvarxjkdiag;
1717 	      SCIP_Bool isauxvarxjldiag;
1718 	      SCIP_Real solxik;
1719 	      SCIP_Real solxil;
1720 	      SCIP_Real solxjk;
1721 	      SCIP_Real solxjl;
1722 	      SCIP_Real det;
1723 	
1724 	      /* get variables of the i-th minor */
1725 	      SCIP_CALL( getMinorVars(sepadata, i, &auxvarxik, &auxvarxil, &auxvarxjk, &auxvarxjl, &isauxvarxikdiag,
1726 	            &isauxvarxildiag, &isauxvarxjkdiag, &isauxvarxjldiag) );
1727 	
1728 	      /* get current solution values */
1729 	      solxik = SCIPvarGetLPSol(auxvarxik);
1730 	      solxil = SCIPvarGetLPSol(auxvarxil);
1731 	      solxjk = SCIPvarGetLPSol(auxvarxjk);
1732 	      solxjl = SCIPvarGetLPSol(auxvarxjl);
1733 	
1734 	      det = solxik * solxjl - solxil * solxjk;
1735 	
1736 	      if( SCIPisFeasZero(scip, det) )
1737 	         continue;
1738 	
1739 	      if( basicvarpos2tableaurow == NULL )
1740 	      {
1741 	         /* allocate memory */
1742 	         SCIP_CALL( SCIPallocBufferArray(scip, &basicvarpos2tableaurow, SCIPgetNLPCols(scip)) );
1743 	         SCIP_CALL( SCIPhashmapCreate(&tableau, SCIPblkmem(scip), SCIPgetNVars(scip)) );
1744 	
1745 	         /* construct basicvar to tableau row map */
1746 	         SCIP_CALL( constructBasicVars2TableauRowMap(scip, basicvarpos2tableaurow) );
1747 	      }
1748 	      assert(tableau != NULL);
1749 	
1750 	      if( SCIPisFeasPositive(scip, det) )
1751 	      {
1752 	         SCIP_CALL( separateDeterminant(scip, sepa, sepadata, auxvarxik, auxvarxil, auxvarxjk, auxvarxjl, &isauxvarxikdiag,
1753 	                  &isauxvarxildiag, &isauxvarxjkdiag, &isauxvarxjldiag, basicvarpos2tableaurow, tableau, result) );
1754 	      }
1755 	      else
1756 	      {
1757 	         assert(SCIPisFeasNegative(scip, det));
1758 	         SCIP_CALL( separateDeterminant(scip, sepa, sepadata, auxvarxil, auxvarxik, auxvarxjl, auxvarxjk, &isauxvarxildiag,
1759 	                  &isauxvarxikdiag, &isauxvarxjldiag, &isauxvarxjkdiag, basicvarpos2tableaurow, tableau, result) );
1760 	      }
1761 	   }
1762 	
1763 	   /* all minors were feasible, so no memory to free */
1764 	   if( basicvarpos2tableaurow == NULL )
1765 	      return SCIP_OKAY;
1766 	
1767 	   /* free memory */
1768 	   for( i = 0; i < SCIPhashmapGetNEntries(tableau); ++i )
1769 	   {
1770 	      SCIP_HASHMAPENTRY* entry;
1771 	
1772 	      entry = SCIPhashmapGetEntry(tableau, i);
1773 	
1774 	      if( entry != NULL )
1775 	      {
1776 	         SCIP_Real* tableaurow;
1777 	
1778 	         tableaurow = (SCIP_Real *) SCIPhashmapEntryGetImage(entry);
1779 	
1780 	         SCIPfreeBufferArrayNull(scip, &tableaurow);
1781 	      }
1782 	   }
1783 	   SCIPhashmapFree(&tableau);
1784 	   SCIPfreeBufferArray(scip, &basicvarpos2tableaurow);
1785 	
1786 	   return SCIP_OKAY;
1787 	}
1788 	
1789 	/*
1790 	 * Callback methods of separator
1791 	 */
1792 	
1793 	/** copy method for separator plugins (called when SCIP copies plugins) */
1794 	static
1795 	SCIP_DECL_SEPACOPY(sepaCopyMinor)
1796 	{  /*lint --e{715}*/
1797 	   assert(scip != NULL);
1798 	   assert(sepa != NULL);
1799 	   assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
1800 	
1801 	   /* call inclusion method of constraint handler */
1802 	   SCIP_CALL( SCIPincludeSepaInterminor(scip) );
1803 	
1804 	   return SCIP_OKAY;
1805 	}
1806 	
1807 	
1808 	/** destructor of separator to free user data (called when SCIP is exiting) */
1809 	static
1810 	SCIP_DECL_SEPAFREE(sepaFreeMinor)
1811 	{  /*lint --e{715}*/
1812 	   SCIP_SEPADATA* sepadata;
1813 	
1814 	   sepadata = SCIPsepaGetData(sepa);
1815 	   assert(sepadata != NULL);
1816 	   assert(sepadata->minors == NULL);
1817 	   assert(sepadata->nminors == 0);
1818 	   assert(sepadata->minorssize == 0);
1819 	
1820 	   /* free separator data */
1821 	   SCIPfreeBlockMemory(scip, &sepadata);
1822 	   SCIPsepaSetData(sepa, NULL);
1823 	
1824 	   return SCIP_OKAY;
1825 	}
1826 	
1827 	
1828 	/** initialization method of separator (called after problem was transformed) */
1829 	static
1830 	SCIP_DECL_SEPAINIT(sepaInitMinor)
1831 	{  /*lint --e{715}*/
1832 	   SCIP_SEPADATA* sepadata;
1833 	
1834 	   /* get separator data */
1835 	   sepadata = SCIPsepaGetData(sepa);
1836 	   assert(sepadata != NULL);
1837 	   assert(sepadata->randnumgen == NULL);
1838 	
1839 	   /* create random number generator */
1840 	   SCIP_CALL( SCIPcreateRandom(scip, &sepadata->randnumgen, DEFAULT_RANDSEED, TRUE) );
1841 	
1842 	   return SCIP_OKAY;
1843 	}
1844 	
1845 	
1846 	/** deinitialization method of separator (called before transformed problem is freed) */
1847 	static
1848 	SCIP_DECL_SEPAEXIT(sepaExitMinor)
1849 	{  /*lint --e{715}*/
1850 	   SCIP_SEPADATA* sepadata;
1851 	
1852 	   /* get separator data */
1853 	   sepadata = SCIPsepaGetData(sepa);
1854 	   assert(sepadata != NULL);
1855 	   assert(sepadata->randnumgen != NULL);
1856 	
1857 	   /* free random number generator */
1858 	   SCIPfreeRandom(scip, &sepadata->randnumgen);
1859 	
1860 	   return SCIP_OKAY;
1861 	}
1862 	
1863 	
1864 	/** solving process initialization method of separator (called when branch and bound process is about to begin) */
1865 	static
1866 	SCIP_DECL_SEPAINITSOL(sepaInitsolMinor)
1867 	{  /*lint --e{715}*/
1868 	   return SCIP_OKAY;
1869 	}
1870 	
1871 	
1872 	/** solving process deinitialization method of separator (called before branch and bound process data is freed) */
1873 	static
1874 	SCIP_DECL_SEPAEXITSOL(sepaExitsolMinor)
1875 	{  /*lint --e{715}*/
1876 	   SCIP_SEPADATA* sepadata;
1877 	
1878 	   sepadata = SCIPsepaGetData(sepa);
1879 	   assert(sepadata != NULL);
1880 	
1881 	   /* clear separation data */
1882 	   SCIP_CALL( sepadataClear(scip, sepadata) );
1883 	
1884 	   return SCIP_OKAY;
1885 	}
1886 	
1887 	
1888 	/** LP solution separation method of separator */
1889 	static
1890 	SCIP_DECL_SEPAEXECLP(sepaExeclpMinor)
1891 	{  /*lint --e{715}*/
1892 	   SCIP_SEPADATA* sepadata;
1893 	   int ncalls;
1894 	   int currentdepth;
1895 	
1896 	   /* need routine to compute eigenvalues/eigenvectors */
1897 	   if( !SCIPisIpoptAvailableIpopt() )
1898 	      return SCIP_OKAY;
1899 	
1900 	   sepadata = SCIPsepaGetData(sepa);
1901 	   assert(sepadata != NULL);
1902 	   currentdepth = SCIPgetDepth(scip);
1903 	   ncalls = SCIPsepaGetNCallsAtNode(sepa);
1904 	
1905 	   /* only call the separator a given number of times at each node */
1906 	   if( (currentdepth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
1907 	      || (currentdepth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
1908 	   {
1909 	      SCIPdebugMsg(scip, "reached round limit for node\n");
1910 	      return SCIP_OKAY;
1911 	   }
1912 	
1913 	   /* try to detect minors */
1914 	   SCIP_CALL( detectMinors(scip, sepadata) );
1915 	
1916 	   /* call separation method */
1917 	   SCIP_CALL( separatePoint(scip, sepa, result) );
1918 	
1919 	   return SCIP_OKAY;
1920 	}
1921 	
1922 	
1923 	/*
1924 	 * separator specific interface methods
1925 	 */
1926 	
1927 	/** creates the minor separator and includes it in SCIP */
1928 	SCIP_RETCODE SCIPincludeSepaInterminor(
1929 	   SCIP*                 scip                /**< SCIP data structure */
1930 	   )
1931 	{
1932 	   SCIP_SEPADATA* sepadata = NULL;
1933 	   SCIP_SEPA* sepa = NULL;
1934 	
1935 	   /* create minor separator data */
1936 	   SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) );
1937 	   BMSclearMemory(sepadata);
1938 	
1939 	   /* include separator */
1940 	   SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
1941 	         SEPA_USESSUBSCIP, SEPA_DELAY,
1942 	         sepaExeclpMinor, NULL,
1943 	         sepadata) );
1944 	
1945 	   assert(sepa != NULL);
1946 	
1947 	   /* set non fundamental callbacks via setter functions */
1948 	   SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyMinor) );
1949 	   SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeMinor) );
1950 	   SCIP_CALL( SCIPsetSepaInit(scip, sepa, sepaInitMinor) );
1951 	   SCIP_CALL( SCIPsetSepaExit(scip, sepa, sepaExitMinor) );
1952 	   SCIP_CALL( SCIPsetSepaInitsol(scip, sepa, sepaInitsolMinor) );
1953 	   SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolMinor) );
1954 	
1955 	   /* add minor separator parameters */
1956 	   SCIP_CALL( SCIPaddBoolParam(scip,
1957 	         "separating/" SEPA_NAME "/usestrengthening",
1958 	         "whether to use strengthened intersection cuts to separate minors",
1959 	         &sepadata->usestrengthening, FALSE, DEFAULT_USESTRENGTHENING, NULL, NULL) );
1960 	
1961 	   SCIP_CALL( SCIPaddBoolParam(scip,
1962 	         "separating/" SEPA_NAME "/usebounds",
1963 	         "whether to also enforce nonegativity bounds of principle minors",
1964 	         &sepadata->usebounds, FALSE, DEFAULT_USEBOUNDS, NULL, NULL) );
1965 	
1966 	   SCIP_CALL( SCIPaddRealParam(scip,
1967 	         "separating/" SEPA_NAME "/mincutviol",
1968 	         "minimum required violation of a cut",
1969 	         &sepadata->mincutviol, FALSE, DEFAULT_MINCUTVIOL, 0.0, SCIP_REAL_MAX, NULL, NULL) );
1970 	
1971 	   SCIP_CALL( SCIPaddIntParam(scip,
1972 	         "separating/" SEPA_NAME "/maxrounds",
1973 	         "maximal number of separation rounds per node (-1: unlimited)",
1974 	         &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
1975 	
1976 	   SCIP_CALL( SCIPaddIntParam(scip,
1977 	         "separating/" SEPA_NAME "/maxroundsroot",
1978 	         "maximal number of separation rounds in the root node (-1: unlimited)",
1979 	         &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
1980 	
1981 	   return SCIP_OKAY;
1982 	}
1983