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_minor.c
26   	 * @ingroup DEFPLUGINS_SEPA
27   	 * @brief  principal minor separator
28   	 * @author Benjamin Mueller
29   	 *
30   	 * @todo detect non-principal minors and use them to derive split cuts
31   	 */
32   	
33   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
34   	
35   	#include <assert.h>
36   	#include <string.h>
37   	
38   	#include "scip/sepa_minor.h"
39   	#include "scip/cons_nonlinear.h"
40   	#include "scip/lapack_calls.h"
41   	
42   	#define SEPA_NAME              "minor"
43   	#define SEPA_DESC              "separator to ensure that 2x2 principal minors of X - xx' are positive semi-definite"
44   	#define SEPA_PRIORITY                 0
45   	#define SEPA_FREQ                    10
46   	#define SEPA_MAXBOUNDDIST           1.0
47   	#define SEPA_USESSUBSCIP          FALSE /**< does the separator use a secondary SCIP instance? */
48   	#define SEPA_DELAY                FALSE /**< should separation method be delayed, if other separators found cuts? */
49   	
50   	#define DEFAULT_MAXMINORSCONST     3000 /**< default constant for the maximum number of minors, i.e., max(const, fac * # quadratic terms) */
51   	#define DEFAULT_MAXMINORSFAC       10.0 /**< default factor for the maximum number of minors, i.e., max(const, fac * # quadratic terms) */
52   	#define DEFAULT_MINCUTVIOL         1e-4 /**< default minimum required violation of a cut */
53   	#define DEFAULT_RANDSEED            157 /**< default random seed */
54   	#define DEFAULT_MAXROUNDS            10 /**< maximal number of separation rounds per node (-1: unlimited) */
55   	#define DEFAULT_MAXROUNDSROOT        -1 /**< maximal number of separation rounds in the root node (-1: unlimited) */
56   	#define DEFAULT_IGNOREPACKINGCONSS TRUE /**< default for ignoring circle packing constraints during minor detection */
57   	
58   	/*
59   	 * Data structures
60   	 */
61   	
62   	/** separator data */
63   	struct SCIP_SepaData
64   	{
65   	   SCIP_VAR**            minors;             /**< variables of 2x2 minors; each minor is stored like (auxvar_x^2,auxvar_y^2,auxvar_xy) */
66   	   int                   nminors;            /**< total number of minors */
67   	   int                   minorssize;         /**< size of minors array */
68   	   int                   maxminorsconst;     /**< constant for the maximum number of minors, i.e., max(const, fac * # quadratic terms) */
69   	   SCIP_Real             maxminorsfac;       /**< factor for the maximum number of minors, i.e., max(const, fac * # quadratic terms) */
70   	   int                   maxrounds;          /**< maximal number of separation rounds per node (-1: unlimited) */
71   	   int                   maxroundsroot;      /**< maximal number of separation rounds in the root node (-1: unlimited) */
72   	   SCIP_Bool             detectedminors;     /**< has minor detection be called? */
73   	   SCIP_Real             mincutviol;         /**< minimum required violation of a cut */
74   	   SCIP_RANDNUMGEN*      randnumgen;         /**< random number generation */
75   	   SCIP_Bool             ignorepackingconss; /**< whether to ignore circle packing constraints during minor detection */
76   	};
77   	
78   	/*
79   	 * Local methods
80   	 */
81   	
82   	/** helper method to store a 2x2 minor in the separation data */
83   	static
84   	SCIP_RETCODE sepadataAddMinor(
85   	   SCIP*                 scip,               /**< SCIP data structure */
86   	   SCIP_SEPADATA*        sepadata,           /**< separator data */
87   	   SCIP_VAR*             x,                  /**< x variable */
88   	   SCIP_VAR*             y,                  /**< y variable */
89   	   SCIP_VAR*             auxvarxx,           /**< auxiliary variable for x*x */
90   	   SCIP_VAR*             auxvaryy,           /**< auxiliary variable for y*y */
91   	   SCIP_VAR*             auxvarxy            /**< auxiliary variable for x*y */
92   	   )
93   	{
94   	   assert(sepadata != NULL);
95   	   assert(x != NULL);
96   	   assert(y != NULL);
97   	   assert(x != y);
98   	   assert(auxvarxx != NULL);
99   	   assert(auxvaryy != NULL);
100  	   assert(auxvarxy != NULL);
101  	   assert(auxvarxx != auxvaryy);
102  	   assert(auxvarxx != auxvarxy);
103  	   assert(auxvaryy != auxvarxy);
104  	
105  	   SCIPdebugMsg(scip, "store 2x2 minor: %s %s %s for x=%s y=%s\n", SCIPvarGetName(auxvarxx), SCIPvarGetName(auxvaryy),
106  	      SCIPvarGetName(auxvarxy), SCIPvarGetName(x), SCIPvarGetName(y));
107  	
108  	   /* reallocate if necessary */
109  	   if( sepadata->minorssize < 5 * (sepadata->nminors + 1) )
110  	   {
111  	      int newsize = SCIPcalcMemGrowSize(scip, 5 * (sepadata->nminors + 1));
112  	      assert(newsize > 5 * (sepadata->nminors + 1));
113  	
114  	      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(sepadata->minors), sepadata->minorssize, newsize) );
115  	      sepadata->minorssize = newsize;
116  	   }
117  	
118  	   /* store minor */
119  	   sepadata->minors[5 * sepadata->nminors] = x;
120  	   sepadata->minors[5 * sepadata->nminors + 1] = y;
121  	   sepadata->minors[5 * sepadata->nminors + 2] = auxvarxx;
122  	   sepadata->minors[5 * sepadata->nminors + 3] = auxvaryy;
123  	   sepadata->minors[5 * sepadata->nminors + 4] = auxvarxy;
124  	   ++(sepadata->nminors);
125  	
126  	   /* capture variables */
127  	   SCIP_CALL( SCIPcaptureVar(scip, x) );
128  	   SCIP_CALL( SCIPcaptureVar(scip, y) );
129  	   SCIP_CALL( SCIPcaptureVar(scip, auxvarxx) );
130  	   SCIP_CALL( SCIPcaptureVar(scip, auxvaryy) );
131  	   SCIP_CALL( SCIPcaptureVar(scip, auxvarxy) );
132  	
133  	   return SCIP_OKAY;
134  	}
135  	
136  	/** helper method to clear separation data */
137  	static
138  	SCIP_RETCODE sepadataClear(
139  	   SCIP*                 scip,               /**< SCIP data structure */
140  	   SCIP_SEPADATA*        sepadata            /**< separator data */
141  	   )
142  	{
143  	   int i;
144  	
145  	   assert(sepadata != NULL);
146  	
147  	   SCIPdebugMsg(scip, "clear separation data\n");
148  	
149  	   /* release captured variables */
150  	   for( i = 0; i < 5 * sepadata->nminors; ++i )
151  	   {
152  	      assert(sepadata->minors[i] != NULL);
153  	      SCIP_CALL( SCIPreleaseVar(scip, &sepadata->minors[i]) );
154  	   }
155  	
156  	   /* free memory */
157  	   SCIPfreeBlockMemoryArrayNull(scip, &sepadata->minors, sepadata->minorssize);
158  	
159  	   /* reset counters */
160  	   sepadata->nminors = 0;
161  	   sepadata->minorssize = 0;
162  	
163  	   return SCIP_OKAY;
164  	}
165  	
166  	/** helper method to identify non-overlapping constraints in circle packing */
167  	static
168  	SCIP_Bool isPackingCons(
169  	   SCIP*                 scip,               /**< SCIP data structure */
170  	   SCIP_CONS*            cons                /**< nonlinear constraint */
171  	   )
172  	{
173  	   SCIP_EXPR* root;
174  	   SCIP_VAR* quadvars[4] = {NULL, NULL, NULL, NULL};
175  	   SCIP_VAR* bilinvars[4] = {NULL, NULL, NULL, NULL};
176  	   int nbilinvars = 0;
177  	   int nquadvars = 0;
178  	   int nchildren;
179  	   int i;
180  	
181  	   assert(scip != NULL);
182  	   assert(cons != NULL);
183  	
184  	   root = SCIPgetExprNonlinear(cons);
185  	   assert(root != NULL);
186  	   nchildren = SCIPexprGetNChildren(root);
187  	
188  	   /* non-overlapping constraint has 6 terms (2 bilinear + 4 quadratic) */
189  	   if( nchildren != 6 || !SCIPisExprSum(scip, root) )
190  	      return FALSE;
191  	
192  	   for( i = 0; i < nchildren; ++i )
193  	   {
194  	      SCIP_EXPR* expr;
195  	      SCIP_EXPR** children;
196  	
197  	      /* get child */
198  	      expr = SCIPexprGetChildren(root)[i];
199  	      assert(expr != NULL);
200  	      children = SCIPexprGetChildren(expr);
201  	
202  	      /* case: expr = x^2; x is no auxiliary variable */
203  	      if( SCIPisExprPower(scip, expr) && SCIPgetExponentExprPow(expr) == 2.0
204  	         && SCIPisExprVar(scip, children[0]) )
205  	      {
206  	         SCIP_VAR* x;
207  	
208  	         /* too many quadratic variables -> stop */
209  	         if( nquadvars > 3 )
210  	            return FALSE;
211  	
212  	         x = SCIPgetVarExprVar(children[0]);
213  	         assert(x != NULL);
214  	
215  	         quadvars[nquadvars++] = x;
216  	      }
217  	      /* case: expr = x * y; x and y are no auxiliary variables */
218  	      else if( SCIPisExprProduct(scip, expr) && SCIPexprGetNChildren(expr) == 2
219  	         && SCIPisExprVar(scip, children[0]) && SCIPisExprVar(scip, children[1]) )
220  	      {
221  	         SCIP_VAR* x;
222  	         SCIP_VAR* y;
223  	
224  	         /* too many bilinear variables -> stop */
225  	         if( nbilinvars > 2 )
226  	            return FALSE;
227  	
228  	         x = SCIPgetVarExprVar(children[0]);
229  	         assert(x != NULL);
230  	         y = SCIPgetVarExprVar(children[1]);
231  	         assert(y != NULL);
232  	         assert(x != y);
233  	
234  	         bilinvars[nbilinvars++] = x;
235  	         bilinvars[nbilinvars++] = y;
236  	      }
237  	      else
238  	      {
239  	         return FALSE;
240  	      }
241  	   }
242  	
243  	   /* number of bilinear and quadratic terms do not fit */
244  	   if( nbilinvars != 4 || nquadvars != 4 )
245  	      return FALSE;
246  	
247  	   /* each quadratic variable has to appear in exactly one bilinear terms */
248  	   for( i = 0; i < nquadvars; ++i )
249  	   {
250  	      int counter = 0;
251  	      int j;
252  	
253  	      for( j = 0; j < nbilinvars; ++j )
254  	      {
255  	         if( quadvars[i] == bilinvars[j] )
256  	            ++counter;
257  	      }
258  	
259  	      if( counter != 1 )
260  	         return FALSE;
261  	   }
262  	
263  	   return TRUE;
264  	}
265  	
266  	/** helper method to get the variables associated to a minor */
267  	static
268  	SCIP_RETCODE getMinorVars(
269  	   SCIP_SEPADATA*        sepadata,           /**< separator data */
270  	   int                   idx,                /**< index of the stored minor */
271  	   SCIP_VAR**            x,                  /**< pointer to store x variable */
272  	   SCIP_VAR**            y,                  /**< pointer to store x variable */
273  	   SCIP_VAR**            auxvarxx,           /**< pointer to store auxiliary variable for x*x */
274  	   SCIP_VAR**            auxvaryy,           /**< pointer to store auxiliary variable for y*y */
275  	   SCIP_VAR**            auxvarxy            /**< pointer to store auxiliary variable for x*y */
276  	   )
277  	{
278  	   assert(sepadata != NULL);
279  	   assert(idx >= 0 && idx < sepadata->nminors);
280  	   assert(auxvarxx != NULL);
281  	   assert(auxvaryy != NULL);
282  	   assert(auxvarxy != NULL);
283  	
284  	   *x = sepadata->minors[5 * idx];
285  	   *y = sepadata->minors[5 * idx + 1];
286  	   *auxvarxx = sepadata->minors[5 * idx + 2];
287  	   *auxvaryy = sepadata->minors[5 * idx + 3];
288  	   *auxvarxy = sepadata->minors[5 * idx + 4];
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* quadmap;
303  	   SCIP_VAR** xs;
304  	   SCIP_VAR** ys;
305  	   SCIP_VAR** auxvars;
306  	   int* perm = NULL;
307  	   int nbilinterms = 0;
308  	   int nquadterms = 0;
309  	   int maxminors;
310  	   int c;
311  	   int i;
312  	
313  	#ifdef SCIP_STATISTIC
314  	   SCIP_Real totaltime = -SCIPgetTotalTime(scip);
315  	#endif
316  	
317  	   assert(sepadata != NULL);
318  	
319  	   /* check whether minor detection has been called already */
320  	   if( sepadata->detectedminors )
321  	      return SCIP_OKAY;
322  	
323  	   assert(sepadata->minors == NULL);
324  	   assert(sepadata->nminors == 0);
325  	
326  	   /* we assume that the auxiliary variables in the nonlinear constraint handler have been already generated */
327  	   sepadata->detectedminors = TRUE;
328  	
329  	   /* check whether there are nonlinear constraints available */
330  	   conshdlr = SCIPfindConshdlr(scip, "nonlinear");
331  	   if( conshdlr == NULL || SCIPconshdlrGetNConss(conshdlr) == 0 )
332  	      return SCIP_OKAY;
333  	
334  	   SCIPdebugMsg(scip, "call detectMinors()\n");
335  	
336  	   /* allocate memory */
337  	   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
338  	   SCIP_CALL( SCIPhashmapCreate(&quadmap, SCIPblkmem(scip), SCIPgetNVars(scip)) );
339  	   SCIP_CALL( SCIPallocBufferArray(scip, &xs, SCIPgetNVars(scip)) );
340  	   SCIP_CALL( SCIPallocBufferArray(scip, &ys, SCIPgetNVars(scip)) );
341  	   SCIP_CALL( SCIPallocBufferArray(scip, &auxvars, SCIPgetNVars(scip)) );
342  	
343  	   /* initialize iterator */
344  	   SCIP_CALL( SCIPexpriterInit(it, NULL, SCIP_EXPRITER_DFS, FALSE) );
345  	   SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_ENTEREXPR);
346  	
347  	   for( c = 0; c < SCIPconshdlrGetNConss(conshdlr); ++c )
348  	   {
349  	      SCIP_CONS* cons;
350  	      SCIP_EXPR* expr;
351  	      SCIP_EXPR* root;
352  	
353  	      cons = SCIPconshdlrGetConss(conshdlr)[c];
354  	      assert(cons != NULL);
355  	      root = SCIPgetExprNonlinear(cons);
356  	      assert(root != NULL);
357  	
358  	      /* ignore circle packing constraints; the motivation for this is that in circle packing instance not only the SDP
359  	       * relaxation is weak (see "Packing circles in a square: a theoretical comparison of various convexification
360  	       * techniques", http://www.optimization-online.org/DB_HTML/2017/03/5911.html), but it also hurts performance
361  	       */
362  	      if( sepadata->ignorepackingconss && isPackingCons(scip, cons) )
363  	      {
364  	         SCIPdebugMsg(scip, "ignore packing constraints %s\n", SCIPconsGetName(cons));
365  	         continue;
366  	      }
367  	
368  	      for( expr = SCIPexpriterRestartDFS(it, root); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) /*lint !e441*/ /*lint !e440*/
369  	      {
370  	         SCIP_EXPR** children;
371  	         SCIP_VAR* auxvar;
372  	
373  	         SCIPdebugMsg(scip, "visit expression %p in constraint %s\n", (void*)expr, SCIPconsGetName(cons));
374  	
375  	         /* check whether the expression has an auxiliary variable */
376  	         auxvar = SCIPgetExprAuxVarNonlinear(expr);
377  	         if( auxvar == NULL )
378  	         {
379  	            SCIPdebugMsg(scip, "expression has no auxiliary variable -> skip\n");
380  	            continue;
381  	         }
382  	
383  	         children = SCIPexprGetChildren(expr);
384  	
385  	         /* check for expr = (x)^2 */
386  	         if( SCIPexprGetNChildren(expr) == 1 && SCIPisExprPower(scip, expr)
387  	            && SCIPgetExponentExprPow(expr) == 2.0
388  	            && SCIPgetExprAuxVarNonlinear(children[0]) != NULL )
389  	         {
390  	            SCIP_VAR* quadvar;
391  	
392  	            assert(children[0] != NULL);
393  	
394  	            quadvar = SCIPgetExprAuxVarNonlinear(children[0]);
395  	            assert(quadvar != NULL);
396  	            assert(!SCIPhashmapExists(quadmap, (void*)quadvar));
397  	            SCIPdebugMsg(scip, "found %s = (%s)^2\n", SCIPvarGetName(auxvar), SCIPvarGetName(quadvar));
398  	
399  	            /* hash the quadratic variable to its corresponding auxiliary variable */
400  	            SCIP_CALL( SCIPhashmapInsert(quadmap, (void*)quadvar, auxvar) );
401  	            ++nquadterms;
402  	         }
403  	         /* check for expr = x * y */
404  	         else if( SCIPexprGetNChildren(expr) == 2 && SCIPisExprProduct(scip, expr)
405  	            && SCIPgetExprAuxVarNonlinear(children[0]) != NULL && SCIPgetExprAuxVarNonlinear(children[1]) != NULL )
406  	         {
407  	            SCIP_VAR* x;
408  	            SCIP_VAR* y;
409  	
410  	            assert(children[0] != NULL);
411  	            assert(children[1] != NULL);
412  	
413  	            x = SCIPgetExprAuxVarNonlinear(children[0]);
414  	            y = SCIPgetExprAuxVarNonlinear(children[1]);
415  	
416  	            /* ignore binary variables */
417  	            if( !SCIPvarIsBinary(x) && !SCIPvarIsBinary(y) )
418  	            {
419  	               xs[nbilinterms] = SCIPgetExprAuxVarNonlinear(children[0]);
420  	               ys[nbilinterms] = SCIPgetExprAuxVarNonlinear(children[1]);
421  	               auxvars[nbilinterms] = auxvar;
422  	               SCIPdebugMsg(scip, "found %s = %s * %s\n", SCIPvarGetName(auxvar), SCIPvarGetName(xs[nbilinterms]), SCIPvarGetName(ys[nbilinterms]));
423  	               ++nbilinterms;
424  	            }
425  	         }
426  	      }
427  	   }
428  	   assert(nbilinterms < SCIPgetNVars(scip));
429  	   SCIPdebugMsg(scip, "stored %d bilinear terms in total\n", nbilinterms);
430  	
431  	   /* use max(maxminorsconst, maxminorsfac * # quadratic terms) as a limit for the maximum number of minors */
432  	   maxminors = (int) MAX(sepadata->maxminorsconst, sepadata->maxminorsfac * nquadterms);
433  	   SCIPdebugMsg(scip, "maximum number of minors = %d\n", maxminors);
434  	
435  	   /* permute bilinear terms if there are too many of them; the motivation for this is that we don't want to
436  	    * prioritize variables because of the order in the bilinear terms where they appear; however, variables that
437  	    * appear more often in bilinear terms might be more important than others so the corresponding bilinear terms
438  	    * are more likely to be chosen
439  	    */
440  	   if( maxminors < nbilinterms && maxminors < SQR(nquadterms) )
441  	   {
442  	      SCIP_CALL( SCIPallocBufferArray(scip, &perm, nbilinterms) );
443  	
444  	      for( i = 0; i < nbilinterms; ++i )
445  	         perm[i] = i;
446  	
447  	      /* permute array */
448  	      SCIPrandomPermuteIntArray(sepadata->randnumgen, perm, 0, nbilinterms);
449  	   }
450  	
451  	   /* store 2x2 principal minors */
452  	   for( i = 0; i < nbilinterms && sepadata->nminors < maxminors; ++i )
453  	   {
454  	      SCIP_VAR* x;
455  	      SCIP_VAR* y;
456  	      SCIP_VAR* auxvarxy;
457  	
458  	      if( perm == NULL )
459  	      {
460  	         x = xs[i];
461  	         y = ys[i];
462  	         auxvarxy = auxvars[i];
463  	      }
464  	      else
465  	      {
466  	         x = xs[perm[i]];
467  	         y = ys[perm[i]];
468  	         auxvarxy = auxvars[perm[i]];
469  	      }
470  	
471  	      assert(x != NULL);
472  	      assert(y != NULL);
473  	      assert(auxvarxy != NULL);
474  	      assert(x != y);
475  	
476  	      if( SCIPhashmapExists(quadmap, (void*)x) && SCIPhashmapExists(quadmap, (void*)y) )
477  	      {
478  	         SCIP_VAR* auxvarxx;
479  	         SCIP_VAR* auxvaryy;
480  	
481  	         auxvarxx = (SCIP_VAR*)SCIPhashmapGetImage(quadmap, (void*)x);
482  	         assert(auxvarxx != NULL);
483  	         auxvaryy = (SCIP_VAR*)SCIPhashmapGetImage(quadmap, (void*)y);
484  	         assert(auxvaryy != NULL);
485  	
486  	         /* store minor into the separation data */
487  	         SCIP_CALL( sepadataAddMinor(scip, sepadata, x, y, auxvarxx, auxvaryy, auxvarxy) );
488  	      }
489  	   }
490  	   SCIPdebugMsg(scip, "found %d principal minors in total\n", sepadata->nminors);
491  	
492  	   /* free memory */
493  	   SCIPfreeBufferArrayNull(scip, &perm);
494  	   SCIPfreeBufferArray(scip, &auxvars);
495  	   SCIPfreeBufferArray(scip, &ys);
496  	   SCIPfreeBufferArray(scip, &xs);
497  	   SCIPhashmapFree(&quadmap);
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  	/** helper method to compute eigenvectors and eigenvalues */
509  	static
510  	SCIP_RETCODE getEigenValues(
511  	   SCIP*                 scip,               /**< SCIP data structure */
512  	   SCIP_Real             x,                  /**< solution value of x */
513  	   SCIP_Real             y,                  /**< solution value of y */
514  	   SCIP_Real             xx,                 /**< solution value of x*x */
515  	   SCIP_Real             yy,                 /**< solution value of y*y */
516  	   SCIP_Real             xy,                 /**< solution value of x*y */
517  	   SCIP_Real*            eigenvals,          /**< array to store eigenvalues (at least of size 3) */
518  	   SCIP_Real*            eigenvecs,          /**< array to store eigenvalues (at least of size 9) */
519  	   SCIP_Bool*            success             /**< pointer to store whether eigenvalue computation was successful */
520  	   )
521  	{
522  	   assert(eigenvals != NULL);
523  	   assert(eigenvecs != NULL);
524  	   assert(success != NULL);
525  	
526  	   *success = TRUE;
527  	
528  	   /* construct matrix */
529  	   eigenvecs[0] = 1.0;
530  	   eigenvecs[1] = x;
531  	   eigenvecs[2] = y;
532  	   eigenvecs[3] = x;
533  	   eigenvecs[4] = xx;
534  	   eigenvecs[5] = xy;
535  	   eigenvecs[6] = y;
536  	   eigenvecs[7] = xy;
537  	   eigenvecs[8] = yy;
538  	
539  	   /* use LAPACK to compute the eigenvalues and eigenvectors */
540  	   if( SCIPlapackComputeEigenvalues(SCIPbuffer(scip), TRUE, 3, eigenvecs, eigenvals) != SCIP_OKAY )
541  	   {
542  	      SCIPdebugMsg(scip, "Failed to compute eigenvalues and eigenvectors of augmented quadratic form matrix.\n");
543  	      *success = FALSE;
544  	   }
545  	
546  	   return SCIP_OKAY;
547  	}
548  	
549  	/** generate and add a cut */
550  	static
551  	SCIP_RETCODE addCut(
552  	   SCIP*                 scip,               /**< SCIP data structure */
553  	   SCIP_SEPA*            sepa,               /**< separator */
554  	   SCIP_SOL*             sol,                /**< solution to separate (might be NULL) */
555  	   SCIP_VAR*             x,                  /**< x variable */
556  	   SCIP_VAR*             y,                  /**< y variable */
557  	   SCIP_VAR*             xx,                 /**< auxiliary variable for x*x */
558  	   SCIP_VAR*             yy,                 /**< auxiliary variable for y*y */
559  	   SCIP_VAR*             xy,                 /**< auxiliary variable for x*y */
560  	   SCIP_Real*            eigenvec,           /**< array containing an eigenvector */
561  	   SCIP_Real             eigenval,           /**< eigenvalue */
562  	   SCIP_Real             mincutviol,         /**< minimal required violation */
563  	   SCIP_RESULT*          result              /**< pointer to update the result */
564  	   )
565  	{
566  	   SCIP_VAR* vars[5] = {x, y, xx, yy, xy};
567  	   SCIP_Real coefs[5];
568  	   SCIP_Real constant;
569  	   SCIP_ROWPREP* rowprep;
570  	   SCIP_Bool success;
571  	
572  	   assert(x != NULL);
573  	   assert(y != NULL);
574  	   assert(xx != NULL);
575  	   assert(yy != NULL);
576  	   assert(xy != NULL);
577  	   assert(eigenvec != NULL);
578  	   assert(mincutviol >= 0.0);
579  	   assert(result != NULL);
580  	
581  	   /* check whether the resulting cut is violated enough */
582  	   if( !SCIPisFeasLT(scip, eigenval, -mincutviol) )
583  	      return SCIP_OKAY;
584  	
585  	   /* the resulting cut reads as
586  	    *              (1 x  y )  (v0)
587  	    *  (v0 v1 v2)  (x xx xy)  (v1)  >= 0
588  	    *              (y xy yy)  (v2)
589  	    *  where v is the eigenvector corresponding to a negative eigenvalue
590  	    *  that is,
591  	    *  v0^2 + 2 v0 v1 * x + 2 v0 v2 * y + v1^2 * xx + v2^2 * yy + 2 v1 v2 * xy >= 0
592  	    */
593  	   constant = SQR(eigenvec[0]);
594  	   coefs[0] = 2.0 * eigenvec[0] * eigenvec[1];
595  	   coefs[1] = 2.0 * eigenvec[0] * eigenvec[2];
596  	   coefs[2] = SQR(eigenvec[1]);
597  	   coefs[3] = SQR(eigenvec[2]);
598  	   coefs[4] = 2.0 * eigenvec[1] * eigenvec[2];
599  	
600  	   /* create rowprep */
601  	   SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, SCIP_SIDETYPE_LEFT, FALSE) );
602  	   SCIP_CALL( SCIPaddRowprepTerms(scip, rowprep, 5, vars, coefs) );
603  	   SCIProwprepAddConstant(rowprep, constant);
604  	   SCIPdebug( SCIPprintRowprep(scip, rowprep, NULL) );
605  	   SCIPdebugMsg(scip, "cut violation %g mincutviol = %g\n", SCIPgetRowprepViolation(scip, rowprep, sol, NULL), mincutviol);
606  	
607  	   /* cleanup coefficient and side, esp treat epsilon to integral values; don't consider scaling up here */
608  	   SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, NULL, 0.0, NULL, &success) );
609  	
610  	   /* check cut violation */
611  	   if( success && SCIPgetRowprepViolation(scip, rowprep, sol, NULL) > mincutviol )
612  	   {
613  	      SCIP_ROW* row;
614  	      SCIP_Bool infeasible;
615  	
616  	      /* set name of rowprep */
617  	      (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "minor_%s_%s_%s_%lld", SCIPvarGetName(xx), SCIPvarGetName(yy),
618  	         SCIPvarGetName(xy), SCIPgetNLPs(scip));
619  	
620  	      /* create, add, and release row */
621  	      SCIP_CALL( SCIPgetRowprepRowSepa(scip, &row, rowprep, sepa) );
622  	      SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
623  	      SCIP_CALL( SCIPreleaseRow(scip, &row) );
624  	
625  	      /* update result pointer */
626  	      *result = infeasible ? SCIP_CUTOFF : SCIP_SEPARATED;
627  	   }
628  	
629  	   /* free rowprep */
630  	   SCIPfreeRowprep(scip, &rowprep);
631  	
632  	   return SCIP_OKAY;
633  	}
634  	
635  	/** separates cuts for stored principal minors */
636  	static
637  	SCIP_RETCODE separatePoint(
638  	   SCIP*                 scip,               /**< SCIP data structure */
639  	   SCIP_SEPA*            sepa,               /**< separator */
640  	   SCIP_SOL*             sol,                /**< primal solution that should be separated, or NULL for LP solution */
641  	   SCIP_RESULT*          result              /**< pointer to store the result of the separation call */
642  	   )
643  	{
644  	   SCIP_SEPADATA* sepadata;
645  	   int i;
646  	
647  	   assert(sepa != NULL);
648  	   assert(result != NULL);
649  	
650  	   *result = SCIP_DIDNOTRUN;
651  	
652  	   sepadata = SCIPsepaGetData(sepa);
653  	   assert(sepadata != NULL);
654  	
655  	   /* check whether there are some minors available */
656  	   if( sepadata->nminors == 0 )
657  	      return SCIP_OKAY;
658  	
659  	   *result = SCIP_DIDNOTFIND;
660  	
661  	   for( i = 0; i < sepadata->nminors && (*result != SCIP_CUTOFF); ++i )
662  	   {
663  	      SCIP_Real eigenvals[3];
664  	      SCIP_Real eigenvecs[9];
665  	      SCIP_VAR* x;
666  	      SCIP_VAR* y;
667  	      SCIP_VAR* xx;
668  	      SCIP_VAR* yy;
669  	      SCIP_VAR* xy;
670  	      SCIP_Real solx;
671  	      SCIP_Real soly;
672  	      SCIP_Real solxx;
673  	      SCIP_Real solyy;
674  	      SCIP_Real solxy;
675  	      SCIP_Bool success;
676  	      int k;
677  	
678  	      /* get variables of the i-th minor */
679  	      SCIP_CALL( getMinorVars(sepadata, i, &x, &y, &xx, &yy, &xy) );
680  	      assert(x != NULL);
681  	      assert(y != NULL);
682  	      assert(xx != NULL);
683  	      assert(yy != NULL);
684  	      assert(xy != NULL);
685  	
686  	      /* get current solution values */
687  	      solx = SCIPgetSolVal(scip, sol, x);
688  	      soly = SCIPgetSolVal(scip, sol, y);
689  	      solxx = SCIPgetSolVal(scip, sol, xx);
690  	      solyy = SCIPgetSolVal(scip, sol, yy);
691  	      solxy = SCIPgetSolVal(scip, sol, xy);
692  	      SCIPdebugMsg(scip, "solution values (x,y,xx,yy,xy)=(%g,%g,%g,%g,%g)\n", solx, soly, solxx, solyy, solxy);
693  	
694  	      /* compute eigenvalues and eigenvectors */
695  	      SCIP_CALL( getEigenValues(scip, solx, soly, solxx, solyy, solxy, eigenvals, eigenvecs, &success) );
696  	      if( !success )
697  	         continue;
698  	
699  	      /* try to generate a cut for each negative eigenvalue */
700  	      for( k = 0; k < 3 && (*result != SCIP_CUTOFF); ++k )
701  	      {
702  	         SCIPdebugMsg(scip, "eigenvalue = %g  eigenvector = (%g,%g,%g)\n", eigenvals[k], eigenvecs[3*k], eigenvecs[3*k + 1], eigenvecs[3*k + 2]);
703  	         SCIP_CALL( addCut(scip, sepa, sol, x, y, xx, yy, xy, &eigenvecs[3*k], eigenvals[k], sepadata->mincutviol, result) );
704  	         SCIPdebugMsg(scip, "result: %d\n", *result);
705  	      }
706  	   }
707  	
708  	   return SCIP_OKAY;
709  	}
710  	
711  	/*
712  	 * Callback methods of separator
713  	 */
714  	
715  	/** copy method for separator plugins (called when SCIP copies plugins) */
716  	static
717  	SCIP_DECL_SEPACOPY(sepaCopyMinor)
718  	{  /*lint --e{715}*/
719  	   assert(scip != NULL);
720  	   assert(sepa != NULL);
721  	   assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
722  	
723  	   /* call inclusion method of constraint handler */
724  	   SCIP_CALL( SCIPincludeSepaMinor(scip) );
725  	
726  	   return SCIP_OKAY;
727  	}
728  	
729  	
730  	/** destructor of separator to free user data (called when SCIP is exiting) */
731  	static
732  	SCIP_DECL_SEPAFREE(sepaFreeMinor)
733  	{  /*lint --e{715}*/
734  	   SCIP_SEPADATA* sepadata;
735  	
736  	   sepadata = SCIPsepaGetData(sepa);
737  	   assert(sepadata != NULL);
738  	   assert(sepadata->minors == NULL);
739  	   assert(sepadata->nminors == 0);
740  	   assert(sepadata->minorssize == 0);
741  	
742  	   /* free separator data */
743  	   SCIPfreeBlockMemory(scip, &sepadata);
744  	   SCIPsepaSetData(sepa, NULL);
745  	
746  	   return SCIP_OKAY;
747  	}
748  	
749  	
750  	/** initialization method of separator (called after problem was transformed) */
751  	static
752  	SCIP_DECL_SEPAINIT(sepaInitMinor)
753  	{  /*lint --e{715}*/
754  	   SCIP_SEPADATA* sepadata;
755  	
756  	   /* get separator data */
757  	   sepadata = SCIPsepaGetData(sepa);
758  	   assert(sepadata != NULL);
759  	   assert(sepadata->randnumgen == NULL);
760  	
761  	   /* create random number generator */
762  	   SCIP_CALL( SCIPcreateRandom(scip, &sepadata->randnumgen, DEFAULT_RANDSEED, TRUE) );
763  	
764  	   return SCIP_OKAY;
765  	}
766  	
767  	
768  	/** deinitialization method of separator (called before transformed problem is freed) */
769  	static
770  	SCIP_DECL_SEPAEXIT(sepaExitMinor)
771  	{  /*lint --e{715}*/
772  	   SCIP_SEPADATA* sepadata;
773  	
774  	   /* get separator data */
775  	   sepadata = SCIPsepaGetData(sepa);
776  	   assert(sepadata != NULL);
777  	   assert(sepadata->randnumgen != NULL);
778  	
779  	   /* free random number generator */
780  	   SCIPfreeRandom(scip, &sepadata->randnumgen);
781  	
782  	   return SCIP_OKAY;
783  	}
784  	
785  	
786  	/** solving process initialization method of separator (called when branch and bound process is about to begin) */
787  	static
788  	SCIP_DECL_SEPAINITSOL(sepaInitsolMinor)
789  	{  /*lint --e{715}*/
790  	   return SCIP_OKAY;
791  	}
792  	
793  	
794  	/** solving process deinitialization method of separator (called before branch and bound process data is freed) */
795  	static
796  	SCIP_DECL_SEPAEXITSOL(sepaExitsolMinor)
797  	{  /*lint --e{715}*/
798  	   SCIP_SEPADATA* sepadata;
799  	
800  	   sepadata = SCIPsepaGetData(sepa);
801  	   assert(sepadata != NULL);
802  	
803  	   /* clear separation data */
804  	   SCIP_CALL( sepadataClear(scip, sepadata) );
805  	
806  	   return SCIP_OKAY;
807  	}
808  	
809  	
810  	/** LP solution separation method of separator */
811  	static
812  	SCIP_DECL_SEPAEXECLP(sepaExeclpMinor)
813  	{  /*lint --e{715}*/
814  	   SCIP_SEPADATA* sepadata;
815  	   int ncalls;
816  	
817  	   /* need routine to compute eigenvalues/eigenvectors */
818  	   if( ! SCIPlapackIsAvailable() )
819  	      return SCIP_OKAY;
820  	
821  	   sepadata = SCIPsepaGetData(sepa);
822  	   assert(sepadata != NULL);
823  	   ncalls = SCIPsepaGetNCallsAtNode(sepa);
824  	
825  	   /* only call the separator a given number of times at each node */
826  	   if( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
827  	      || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
828  	   {
829  	      SCIPdebugMsg(scip, "reached round limit for node\n");
830  	      return SCIP_OKAY;
831  	   }
832  	
833  	   /* try to detect minors */
834  	   SCIP_CALL( detectMinors(scip, sepadata) );
835  	
836  	   /* call separation method */
837  	   SCIP_CALL( separatePoint(scip, sepa, NULL, result) );
838  	
839  	   return SCIP_OKAY;
840  	}
841  	
842  	
843  	/** arbitrary primal solution separation method of separator */
844  	static
845  	SCIP_DECL_SEPAEXECSOL(sepaExecsolMinor)
846  	{  /*lint --e{715}*/
847  	   SCIP_SEPADATA* sepadata;
848  	   int ncalls;
849  	
850  	   /* need routine to compute eigenvalues/eigenvectors */
851  	   if( ! SCIPlapackIsAvailable() )
852  	      return SCIP_OKAY;
853  	
854  	   sepadata = SCIPsepaGetData(sepa);
855  	   assert(sepadata != NULL);
856  	   ncalls = SCIPsepaGetNCallsAtNode(sepa);
857  	
858  	   /* only call the separator a given number of times at each node */
859  	   if( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
860  	      || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
861  	   {
862  	      SCIPdebugMsg(scip, "reached round limit for node\n");
863  	      return SCIP_OKAY;
864  	   }
865  	
866  	   /* try to detect minors */
867  	   SCIP_CALL( detectMinors(scip, SCIPsepaGetData(sepa)) );
868  	
869  	   /* call separation method */
870  	   SCIP_CALL( separatePoint(scip, sepa, sol, result) );
871  	
872  	   return SCIP_OKAY;
873  	}
874  	
875  	/*
876  	 * separator specific interface methods
877  	 */
878  	
879  	/** creates the minor separator and includes it in SCIP */
880  	SCIP_RETCODE SCIPincludeSepaMinor(
881  	   SCIP*                 scip                /**< SCIP data structure */
882  	   )
883  	{
884  	   SCIP_SEPADATA* sepadata = NULL;
885  	   SCIP_SEPA* sepa = NULL;
886  	
887  	   /* create minor separator data */
888  	   SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) );
889  	   BMSclearMemory(sepadata);
890  	
891  	   /* include separator */
892  	   SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
893  	         SEPA_USESSUBSCIP, SEPA_DELAY,
894  	         sepaExeclpMinor, sepaExecsolMinor,
895  	         sepadata) );
896  	
897  	   assert(sepa != NULL);
898  	
899  	   /* set non fundamental callbacks via setter functions */
900  	   SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyMinor) );
901  	   SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeMinor) );
902  	   SCIP_CALL( SCIPsetSepaInit(scip, sepa, sepaInitMinor) );
903  	   SCIP_CALL( SCIPsetSepaExit(scip, sepa, sepaExitMinor) );
904  	   SCIP_CALL( SCIPsetSepaInitsol(scip, sepa, sepaInitsolMinor) );
905  	   SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolMinor) );
906  	
907  	   /* add minor separator parameters */
908  	   SCIP_CALL( SCIPaddIntParam(scip,
909  	         "separating/" SEPA_NAME "/maxminorsconst",
910  	         "constant for the maximum number of minors, i.e., max(const, fac * # quadratic terms)",
911  	         &sepadata->maxminorsconst, FALSE, DEFAULT_MAXMINORSCONST, 0, INT_MAX, NULL, NULL) );
912  	
913  	   SCIP_CALL( SCIPaddRealParam(scip,
914  	         "separating/" SEPA_NAME "/maxminorsfac",
915  	         "factor for the maximum number of minors, i.e., max(const, fac * # quadratic terms)",
916  	         &sepadata->maxminorsfac, FALSE, DEFAULT_MAXMINORSFAC, 0.0, SCIP_REAL_MAX, NULL, NULL) );
917  	
918  	   SCIP_CALL( SCIPaddRealParam(scip,
919  	         "separating/" SEPA_NAME "/mincutviol",
920  	         "minimum required violation of a cut",
921  	         &sepadata->mincutviol, FALSE, DEFAULT_MINCUTVIOL, 0.0, SCIP_REAL_MAX, NULL, NULL) );
922  	
923  	   SCIP_CALL( SCIPaddIntParam(scip,
924  	         "separating/" SEPA_NAME "/maxrounds",
925  	         "maximal number of separation rounds per node (-1: unlimited)",
926  	         &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
927  	
928  	   SCIP_CALL( SCIPaddIntParam(scip,
929  	         "separating/" SEPA_NAME "/maxroundsroot",
930  	         "maximal number of separation rounds in the root node (-1: unlimited)",
931  	         &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
932  	
933  	   SCIP_CALL( SCIPaddBoolParam(scip,
934  	         "separating/" SEPA_NAME "/ignorepackingconss",
935  	         "whether to ignore circle packing constraints during minor detection",
936  	         &sepadata->ignorepackingconss, FALSE, DEFAULT_IGNOREPACKINGCONSS, NULL, NULL) );
937  	
938  	   return SCIP_OKAY;
939  	}
940