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_eccuts.c
26   	 * @ingroup DEFPLUGINS_SEPA
27   	 * @brief  edge concave cut separator
28   	 * @author Benjamin Mueller
29   	 */
30   	
31   	/**@todo only count number of fixed variables in the edge concave terms */
32   	/**@todo only add nonlinear row aggregations where at least ...% of the variables (bilinear terms?) are in edge concave
33   	 * terms */
34   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
35   	
36   	#include <assert.h>
37   	#include <string.h>
38   	
39   	#include "scip/scipdefplugins.h"
40   	#include "scip/sepa_eccuts.h"
41   	#include "scip/cons_xor.h"
42   	#include "scip/nlp.h"
43   	#include "tclique/tclique.h"
44   	
45   	#define SEPA_NAME                            "eccuts"
46   	#define SEPA_DESC                            "separator for edge-concave functions"
47   	#define SEPA_PRIORITY            -13000
48   	#define SEPA_FREQ                    -1
49   	#define SEPA_MAXBOUNDDIST           1.0
50   	#define SEPA_USESSUBSCIP          FALSE /**< does the separator use a secondary SCIP instance? */
51   	#define SEPA_DELAY                FALSE /**< should separation method be delayed, if other separators found cuts? */
52   	
53   	#define CLIQUE_MAXFIRSTNODEWEIGHT  1000 /**< maximum weight of branching nodes in level 0; 0 if not used for cliques
54   	                                         *   with at least one fractional node) */
55   	#define CLIQUE_MINWEIGHT              0 /**< lower bound for weight of generated cliques */
56   	#define CLIQUE_MAXNTREENODES      10000 /**< maximal number of nodes of b&b tree */
57   	#define CLIQUE_BACKTRACKFREQ      10000 /**< frequency to backtrack to first level of tree (0: no premature backtracking) */
58   	
59   	#define DEFAULT_DYNAMICCUTS        TRUE /**< should generated cuts be removed from the LP if they are no longer tight? */
60   	#define DEFAULT_MAXROUNDS            10 /**< maximal number of separation rounds per node (-1: unlimited) */
61   	#define DEFAULT_MAXROUNDSROOT       250 /**< maximal number of separation rounds in the root node (-1: unlimited) */
62   	#define DEFAULT_MAXDEPTH             -1 /**< maximal depth at which the separator is applied */
63   	#define DEFAULT_MAXSEPACUTS          10 /**< maximal number of e.c. cuts separated per separation round */
64   	#define DEFAULT_MAXSEPACUTSROOT      50 /**< maximal number of e.c. cuts separated per separation round in root node */
65   	#define DEFAULT_CUTMAXRANGE        1e+7 /**< maximal coefficient range of a cut (maximal coefficient divided by minimal
66   	                                         *   coefficient) in order to be added to LP relaxation */
67   	#define DEFAULT_MINVIOLATION        0.3 /**< minimal violation of an e.c. cut to be separated */
68   	#define DEFAULT_MINAGGRSIZE           3 /**< search for e.c. aggregation of at least this size (has to be >= 3) */
69   	#define DEFAULT_MAXAGGRSIZE           4 /**< search for e.c. aggregation of at most this size (has to be >= minaggrsize) */
70   	#define DEFAULT_MAXBILINTERMS       500 /**< maximum number of bilinear terms allowed to be in a quadratic constraint */
71   	#define DEFAULT_MAXSTALLROUNDS        5 /**< maximum number of unsuccessful rounds in the e.c. aggregation search */
72   	
73   	#define SUBSCIP_NODELIMIT         100LL /**< node limit to solve the sub-SCIP */
74   	
75   	#define ADJUSTFACETTOL             1e-6 /**< adjust resulting facets in checkRikun() up to a violation of this value */
76   	#define USEDUALSIMPLEX             TRUE /**< use dual or primal simplex algorithm? */
77   	
78   	/** first values for \f$2^n\f$ */
79   	static const int poweroftwo[] = { 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192 };
80   	
81   	/*
82   	 * Data structures
83   	 */
84   	
85   	/** data to store a single edge-concave aggregations; an edge-concave aggregation of a quadratic constraint is a subset
86   	 *  of nonconvex bilinear terms
87   	 */
88   	struct EcAggr
89   	{
90   	   SCIP_VAR**            vars;               /**< variables */
91   	   int                   nvars;              /**< number of variables */
92   	   int                   varsize;            /**< size of vars array */
93   	
94   	   SCIP_Real*            termcoefs;          /**< coefficients of bilinear terms */
95   	   int*                  termvars1;          /**< index of the first variable of each bilinear term */
96   	   int*                  termvars2;          /**< index of the second variable of each bilinear term*/
97   	   int                   nterms;             /**< number of bilinear terms in the aggregation */
98   	   int                   termsize;           /**< size of term{coefs,vars1,vars2} arrays */
99   	};
100  	typedef struct EcAggr SCIP_ECAGGR;
101  	
102  	/** data to store all edge-concave aggregations and the remaining part of a nonlinear row of the form g(x) <= rhs */
103  	struct NlrowAggr
104  	{
105  	   SCIP_NLROW*           nlrow;              /**< nonlinear row aggregation */
106  	   SCIP_Bool             rhsaggr;            /**< consider nonlinear row aggregation for g(x) <= rhs (TRUE) or
107  	                                              *   g(x) >= lhs (FALSE) */
108  	
109  	   SCIP_ECAGGR**         ecaggr;             /**< array with all edge-concave aggregations */
110  	   int                   necaggr;            /**< number of edge-concave aggregation */
111  	
112  	   SCIP_VAR**            linvars;            /**< linear variables */
113  	   SCIP_Real*            lincoefs;           /**< linear coefficients */
114  	   int                   nlinvars;           /**< number of linear variables */
115  	   int                   linvarssize;        /**< size of linvars array */
116  	
117  	   SCIP_VAR**            quadvars;           /**< quadratic variables */
118  	   int*                  quadvar2aggr;       /**< stores in which edge-concave aggregation the i-th quadratic variable
119  	                                              *   is contained (< 0: in no edge-concave aggregation) */
120  	   int                   nquadvars;          /**< number of quadratic variables */
121  	   int                   quadvarssize;       /**< size of quadvars array */
122  	
123  	   SCIP_VAR**            remtermvars1;       /**< first quadratic variable of remaining bilinear terms */
124  	   SCIP_VAR**            remtermvars2;       /**< second quadratic variable of remaining bilinear terms */
125  	   SCIP_Real*            remtermcoefs;       /**< coefficients for each remaining bilinear term */
126  	   int                   nremterms;          /**< number of remaining bilinear terms */
127  	   int                   remtermsize;        /**< size of remterm* arrays */
128  	
129  	   SCIP_Real             rhs;                /**< rhs of the nonlinear row */
130  	   SCIP_Real             constant;           /**< constant part of the nonlinear row */
131  	};
132  	typedef struct NlrowAggr SCIP_NLROWAGGR;
133  	
134  	/** separator data */
135  	struct SCIP_SepaData
136  	{
137  	   SCIP_NLROWAGGR**      nlrowaggrs;         /**< array containing all nonlinear row aggregations */
138  	   int                   nnlrowaggrs;        /**< number of nonlinear row aggregations */
139  	   int                   nlrowaggrssize;     /**< size of nlrowaggrs array */
140  	   SCIP_Bool             searchedforaggr;    /**< flag if we already searched for nlrow aggregation candidates */
141  	   int                   minaggrsize;        /**< only search for e.c. aggregations of at least this size (has to be >= 3) */
142  	   int                   maxaggrsize;        /**< only search for e.c. aggregations of at most this size (has to be >= minaggrsize) */
143  	   int                   maxecsize;          /**< largest edge concave aggregation size */
144  	   int                   maxbilinterms;      /**< maximum number of bilinear terms allowed to be in a quadratic constraint */
145  	   int                   maxstallrounds;     /**< maximum number of unsuccessful rounds in the e.c. aggregation search */
146  	
147  	   SCIP_LPI*             lpi;                /**< LP interface to solve the LPs to compute the facets of the convex envelopes */
148  	   int                   lpisize;            /**< maximum size of e.c. aggregations which can be handled by the LP interface */
149  	
150  	   SCIP_Real             cutmaxrange;        /**< maximal coef range of a cut (maximal coefficient divided by minimal
151  	                                              *   coefficient) in order to be added to LP relaxation */
152  	   SCIP_Bool             dynamiccuts;        /**< should generated cuts be removed from the LP if they are no longer tight? */
153  	   SCIP_Real             minviolation;       /**< minimal violation of an e.c. cut to be separated */
154  	
155  	   int                   maxrounds;          /**< maximal number of separation rounds per node (-1: unlimited) */
156  	   int                   maxroundsroot;      /**< maximal number of separation rounds in the root node (-1: unlimited) */
157  	   int                   maxdepth;           /**< maximal depth at which the separator is applied */
158  	   int                   maxsepacuts;        /**< maximal number of e.c. cuts separated per separation round */
159  	   int                   maxsepacutsroot;    /**< maximal number of e.c. cuts separated per separation round in root node */
160  	
161  	#ifdef SCIP_STATISTIC
162  	   SCIP_Real             aggrsearchtime;     /**< total time spent for searching edge concave aggregations */
163  	   int                   nlhsnlrowaggrs;     /**< number of found nonlinear row aggregations for SCIP_NLROWs of the form g(x) <= rhs */
164  	   int                   nrhsnlrowaggrs;     /**< number of found nonlinear row aggregations for SCIP_NLROWs of the form g(x) >= lhs */
165  	#endif
166  	};
167  	
168  	
169  	/*
170  	 * Local methods
171  	 */
172  	
173  	/** creates an empty edge-concave aggregation (without bilinear terms) */
174  	static
175  	SCIP_RETCODE ecaggrCreateEmpty(
176  	   SCIP*                 scip,               /**< SCIP data structure */
177  	   SCIP_ECAGGR**         ecaggr,             /**< pointer to store the edge-concave aggregation */
178  	   int                   nquadvars,          /**< number of quadratic variables */
179  	   int                   nquadterms          /**< number of bilinear terms */
180  	   )
181  	{
182  	   assert(scip != NULL);
183  	   assert(ecaggr != NULL);
184  	   assert(nquadvars > 0);
185  	   assert(nquadterms >= nquadvars);
186  	
187  	   SCIP_CALL( SCIPallocBlockMemory(scip, ecaggr) );
188  	
189  	   (*ecaggr)->nvars = 0;
190  	   (*ecaggr)->nterms = 0;
191  	   (*ecaggr)->varsize = nquadvars;
192  	   (*ecaggr)->termsize = nquadterms;
193  	
194  	   /* allocate enough memory for the quadratic variables and bilinear terms */
195  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*ecaggr)->vars, nquadvars) );
196  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*ecaggr)->termcoefs, nquadterms) );
197  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*ecaggr)->termvars1, nquadterms) );
198  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*ecaggr)->termvars2, nquadterms) );
199  	
200  	   return SCIP_OKAY;
201  	}
202  	
203  	/** frees an edge-concave aggregation */
204  	static
205  	SCIP_RETCODE ecaggrFree(
206  	   SCIP*                 scip,               /**< SCIP data structure */
207  	   SCIP_ECAGGR**         ecaggr              /**< pointer to store the edge-concave aggregation */
208  	   )
209  	{
210  	   assert(scip != NULL);
211  	   assert(ecaggr != NULL);
212  	
213  	   SCIPfreeBlockMemoryArray(scip, &((*ecaggr)->termcoefs), (*ecaggr)->termsize);
214  	   SCIPfreeBlockMemoryArray(scip, &((*ecaggr)->termvars1), (*ecaggr)->termsize);
215  	   SCIPfreeBlockMemoryArray(scip, &((*ecaggr)->termvars2), (*ecaggr)->termsize);
216  	   SCIPfreeBlockMemoryArray(scip, &((*ecaggr)->vars), (*ecaggr)->varsize);
217  	
218  	   SCIPfreeBlockMemory(scip, ecaggr);
219  	   *ecaggr = NULL;
220  	
221  	   return SCIP_OKAY;
222  	}
223  	
224  	/** adds a quadratic variable to an edge-concave aggregation */
225  	static
226  	SCIP_RETCODE ecaggrAddQuadvar(
227  	   SCIP_ECAGGR*          ecaggr,             /**< pointer to store the edge-concave aggregation */
228  	   SCIP_VAR*             x                   /**< first variable */
229  	   )
230  	{
231  	   ecaggr->vars[ ecaggr->nvars++ ] = x;
232  	   return SCIP_OKAY;
233  	}
234  	
235  	/** adds a bilinear term to an edge-concave aggregation */
236  	static
237  	SCIP_RETCODE ecaggrAddBilinTerm(
238  	   SCIP*                 scip,               /**< SCIP data structure */
239  	   SCIP_ECAGGR*          ecaggr,             /**< pointer to store the edge-concave aggregation */
240  	   SCIP_VAR*             x,                  /**< first variable */
241  	   SCIP_VAR*             y,                  /**< second variable */
242  	   SCIP_Real             coef                /**< bilinear coefficient */
243  	   )
244  	{
245  	   int idx1;
246  	   int idx2;
247  	   int i;
248  	
249  	   assert(x != NULL);
250  	   assert(y != NULL);
251  	   assert(ecaggr->nterms + 1 <= ((ecaggr->nvars + 1) * ecaggr->nvars) / 2);
252  	   assert(!SCIPisZero(scip, coef));
253  	
254  	   idx1 = -1;
255  	   idx2 = -1;
256  	
257  	   /* search for the quadratic variables in the e.c. aggregation */
258  	   for( i = 0; i < ecaggr->nvars && (idx1 == -1 || idx2 == -1); ++i )
259  	   {
260  	      if( ecaggr->vars[i] == x )
261  	         idx1 = i;
262  	      if( ecaggr->vars[i] == y )
263  	         idx2 = i;
264  	   }
265  	
266  	   assert(idx1 != -1 && idx2 != -1);
267  	
268  	   ecaggr->termcoefs[ ecaggr->nterms ] = coef;
269  	   ecaggr->termvars1[ ecaggr->nterms ] = idx1;
270  	   ecaggr->termvars2[ ecaggr->nterms ] = idx2;
271  	   ++(ecaggr->nterms);
272  	
273  	   return SCIP_OKAY;
274  	}
275  	
276  	#ifdef SCIP_DEBUG
277  	/** prints an edge-concave aggregation */
278  	static
279  	void ecaggrPrint(
280  	   SCIP*                 scip,               /**< SCIP data structure */
281  	   SCIP_ECAGGR*          ecaggr              /**< pointer to store the edge-concave aggregation */
282  	   )
283  	{
284  	   int i;
285  	
286  	   assert(scip != NULL);
287  	   assert(ecaggr != NULL);
288  	
289  	   SCIPdebugMsg(scip, " nvars = %d nterms = %d\n", ecaggr->nvars, ecaggr->nterms);
290  	   SCIPdebugMsg(scip, " vars: ");
291  	   for( i = 0; i < ecaggr->nvars; ++i )
292  	      SCIPdebugMsgPrint(scip, "%s ", SCIPvarGetName(ecaggr->vars[i]));
293  	   SCIPdebugMsgPrint(scip, "\n");
294  	
295  	   SCIPdebugMsg(scip, " terms: ");
296  	   for( i = 0; i < ecaggr->nterms; ++i )
297  	   {
298  	      SCIP_VAR* x;
299  	      SCIP_VAR* y;
300  	
301  	      x = ecaggr->vars[ ecaggr->termvars1[i] ];
302  	      y = ecaggr->vars[ ecaggr->termvars2[i] ];
303  	      SCIPdebugMsgPrint(scip, "%e %s * %s  ", ecaggr->termcoefs[i], SCIPvarGetName(x), SCIPvarGetName(y) );
304  	   }
305  	   SCIPdebugMsgPrint(scip, "\n");
306  	}
307  	#endif
308  	
309  	/** stores linear terms in a given nonlinear row aggregation */
310  	static
311  	SCIP_RETCODE nlrowaggrStoreLinearTerms(
312  	   SCIP*                 scip,               /**< SCIP data structure */
313  	   SCIP_NLROWAGGR*       nlrowaggr,          /**< nonlinear row aggregation */
314  	   SCIP_VAR**            linvars,            /**< linear variables */
315  	   SCIP_Real*            lincoefs,           /**< linear coefficients */
316  	   int                   nlinvars            /**< number of linear variables */
317  	   )
318  	{
319  	   assert(scip != NULL);
320  	   assert(nlrowaggr != NULL);
321  	   assert(linvars != NULL || nlinvars == 0);
322  	   assert(lincoefs != NULL || nlinvars == 0);
323  	   assert(nlinvars >= 0);
324  	
325  	   nlrowaggr->nlinvars = 0;
326  	   nlrowaggr->linvarssize = 0;
327  	   nlrowaggr->linvars = NULL;
328  	   nlrowaggr->lincoefs = NULL;
329  	
330  	   if( nlinvars == 0 )
331  	      return SCIP_OKAY;
332  	
333  	   SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &nlrowaggr->linvars, linvars, nlinvars) );
334  	   SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &nlrowaggr->lincoefs, lincoefs, nlinvars) );
335  	   nlrowaggr->nlinvars = nlinvars;
336  	   nlrowaggr->linvarssize = nlinvars;
337  	
338  	   /* if we have a nlrow of the form g(x) >= lhs, multiply every coefficient by -1 */
339  	   if( !nlrowaggr->rhsaggr )
340  	   {
341  	      int i;
342  	
343  	      for( i = 0; i < nlrowaggr->nlinvars; ++i )
344  	         nlrowaggr->lincoefs[i] *= -1.0;
345  	   }
346  	
347  	   return SCIP_OKAY;
348  	}
349  	
350  	/** adds linear term to a given nonlinear row aggregation */
351  	static
352  	SCIP_RETCODE nlrowaggrAddLinearTerm(
353  	   SCIP*                 scip,               /**< SCIP data structure */
354  	   SCIP_NLROWAGGR*       nlrowaggr,          /**< nonlinear row aggregation */
355  	   SCIP_VAR*             linvar,             /**< linear variable */
356  	   SCIP_Real             lincoef             /**< coefficient */
357  	   )
358  	{
359  	   assert(scip != NULL);
360  	   assert(nlrowaggr != NULL);
361  	   assert(linvar != NULL);
362  	
363  	   if( nlrowaggr->nlinvars == nlrowaggr->linvarssize )
364  	   {
365  	      int newsize = SCIPcalcMemGrowSize(scip, nlrowaggr->linvarssize+1);
366  	      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &nlrowaggr->linvars, nlrowaggr->linvarssize, newsize) );
367  	      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &nlrowaggr->lincoefs, nlrowaggr->linvarssize, newsize) );
368  	      nlrowaggr->linvarssize = newsize;
369  	   }
370  	   assert(nlrowaggr->linvarssize > nlrowaggr->nlinvars);
371  	
372  	   /* if we have a nlrow of the form g(x) >= lhs, multiply coefficient by -1 */
373  	   if( !nlrowaggr->rhsaggr )
374  	      lincoef = -lincoef;
375  	
376  	   nlrowaggr->linvars[nlrowaggr->nlinvars] = linvar;
377  	   nlrowaggr->lincoefs[nlrowaggr->nlinvars] = lincoef;
378  	   ++nlrowaggr->nlinvars;
379  	
380  	   return SCIP_OKAY;
381  	}
382  	
383  	/** adds quadratic variable to a given nonlinear row aggregation */
384  	static
385  	SCIP_RETCODE nlrowaggrAddQuadraticVar(
386  	   SCIP*                 scip,               /**< SCIP data structure */
387  	   SCIP_NLROWAGGR*       nlrowaggr,          /**< nonlinear row aggregation */
388  	   SCIP_VAR*             quadvar             /**< quadratic variable */
389  	   )
390  	{
391  	   assert(scip != NULL);
392  	   assert(nlrowaggr != NULL);
393  	   assert(quadvar != NULL);
394  	
395  	   SCIP_CALL( SCIPensureBlockMemoryArray(scip, &nlrowaggr->quadvars, &nlrowaggr->quadvarssize, nlrowaggr->nquadvars+1) );
396  	   assert(nlrowaggr->quadvarssize > nlrowaggr->nquadvars);
397  	   nlrowaggr->quadvars[nlrowaggr->nquadvars] = quadvar;
398  	   ++nlrowaggr->nquadvars;
399  	
400  	   return SCIP_OKAY;
401  	}
402  	
403  	/** adds a remaining bilinear term to a given nonlinear row aggregation */
404  	static
405  	SCIP_RETCODE nlrowaggrAddRemBilinTerm(
406  	   SCIP_NLROWAGGR*       nlrowaggr,          /**< nonlinear row aggregation */
407  	   SCIP_VAR*             x,                  /**< first variable */
408  	   SCIP_VAR*             y,                  /**< second variable */
409  	   SCIP_Real             coef                /**< bilinear coefficient */
410  	   )
411  	{
412  	   assert(nlrowaggr != NULL);
413  	   assert(x != NULL);
414  	   assert(y != NULL);
415  	   assert(coef != 0.0);
416  	   assert(nlrowaggr->remtermcoefs != NULL);
417  	   assert(nlrowaggr->remtermvars1 != NULL);
418  	   assert(nlrowaggr->remtermvars2 != NULL);
419  	
420  	   nlrowaggr->remtermcoefs[ nlrowaggr->nremterms ] = coef;
421  	   nlrowaggr->remtermvars1[ nlrowaggr->nremterms ] = x;
422  	   nlrowaggr->remtermvars2[ nlrowaggr->nremterms ] = y;
423  	   ++(nlrowaggr->nremterms);
424  	
425  	   return SCIP_OKAY;
426  	}
427  	
428  	/** creates a nonlinear row aggregation */
429  	static
430  	SCIP_RETCODE nlrowaggrCreate(
431  	   SCIP*                 scip,               /**< SCIP data structure */
432  	   SCIP_NLROW*           nlrow,              /**< nonlinear row */
433  	   SCIP_NLROWAGGR**      nlrowaggr,          /**< pointer to store the nonlinear row aggregation */
434  	   int*                  quadvar2aggr,       /**< mapping between quadratic variables and edge-concave aggregation
435  	                                              *   stores a negative value if the quadratic variables does not belong
436  	                                              *   to any aggregation */
437  	   int                   nfound,             /**< number of edge-concave aggregations */
438  	   SCIP_Bool             rhsaggr             /**< consider nonlinear row aggregation for g(x) <= rhs (TRUE) or
439  	                                              *   lhs <= g(x) (FALSE) */
440  	   )
441  	{
442  	   SCIP_EXPR* expr;
443  	   int* aggrnvars;  /* count the number of variables in each e.c. aggregations */
444  	   int* aggrnterms; /* count the number of bilinear terms in each e.c. aggregations */
445  	   int nquadvars;
446  	   int nremterms;
447  	   int i;
448  	
449  	   assert(scip != NULL);
450  	   assert(nlrow != NULL);
451  	   assert(nlrowaggr != NULL);
452  	   assert(quadvar2aggr != NULL);
453  	   assert(nfound > 0);
454  	
455  	   expr = SCIPnlrowGetExpr(nlrow);
456  	   SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadvars, NULL, NULL, NULL);
457  	   nremterms = 0;
458  	
459  	   SCIP_CALL( SCIPallocClearBufferArray(scip, &aggrnvars, nfound) );
460  	   SCIP_CALL( SCIPallocClearBufferArray(scip, &aggrnterms, nfound) );
461  	
462  	   /* create an empty nonlinear row aggregation */
463  	   SCIP_CALL( SCIPallocBlockMemory(scip, nlrowaggr) );
464  	   (*nlrowaggr)->nlrow = nlrow;
465  	   (*nlrowaggr)->rhsaggr = rhsaggr;
466  	   (*nlrowaggr)->rhs = rhsaggr ? SCIPnlrowGetRhs(nlrow) : -SCIPnlrowGetLhs(nlrow);
467  	   (*nlrowaggr)->constant = rhsaggr ? SCIPnlrowGetConstant(nlrow) : -SCIPnlrowGetConstant(nlrow);
468  	
469  	   (*nlrowaggr)->quadvars = NULL;
470  	   (*nlrowaggr)->nquadvars = 0;
471  	   (*nlrowaggr)->quadvarssize = 0;
472  	   (*nlrowaggr)->quadvar2aggr = NULL;
473  	   (*nlrowaggr)->remtermcoefs = NULL;
474  	   (*nlrowaggr)->remtermvars1 = NULL;
475  	   (*nlrowaggr)->remtermvars2 = NULL;
476  	   (*nlrowaggr)->nremterms = 0;
477  	
478  	   /* copy quadvar2aggr array */
479  	   SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlrowaggr)->quadvar2aggr, quadvar2aggr, nquadvars) );
480  	
481  	   /* store all linear terms */
482  	   SCIP_CALL( nlrowaggrStoreLinearTerms(scip, *nlrowaggr, SCIPnlrowGetLinearVars(nlrow), SCIPnlrowGetLinearCoefs(nlrow),
483  	         SCIPnlrowGetNLinearVars(nlrow)) );
484  	
485  	   /* store all quadratic variables and additional linear terms */
486  	   /* count the number of variables in each e.c. aggregation */
487  	   /* count the number of square and bilinear terms in each e.c. aggregation */
488  	   for( i = 0; i < nquadvars; ++i )
489  	   {
490  	      SCIP_EXPR* qterm;
491  	      SCIP_Real lincoef;
492  	      SCIP_Real sqrcoef;
493  	      int idx1;
494  	      int nadjbilin;
495  	      int* adjbilin;
496  	      int j;
497  	
498  	      SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, &lincoef, &sqrcoef, &nadjbilin, &adjbilin, NULL);
499  	      assert(SCIPisExprVar(scip, qterm));
500  	
501  	      SCIP_CALL( nlrowaggrAddQuadraticVar(scip, *nlrowaggr, SCIPgetVarExprVar(qterm)) );
502  	
503  	      if( lincoef != 0.0 )
504  	      {
505  	         SCIP_CALL( nlrowaggrAddLinearTerm(scip, *nlrowaggr, SCIPgetVarExprVar(qterm), lincoef) );
506  	      }
507  	
508  	      if( quadvar2aggr[i] >= 0)
509  	         ++aggrnvars[ quadvar2aggr[i] ];
510  	
511  	      idx1 = quadvar2aggr[i];
512  	      if( rhsaggr )
513  	         sqrcoef = -sqrcoef;
514  	
515  	      /* variable has to belong to an e.c. aggregation; square term has to be concave */
516  	      if( idx1 >= 0 && SCIPisNegative(scip, sqrcoef) )
517  	         ++aggrnterms[idx1];
518  	      else
519  	         ++nremterms;
520  	
521  	      for( j = 0; j < nadjbilin; ++j )
522  	      {
523  	         SCIP_EXPR* qterm1;
524  	         int pos2;
525  	         int idx2;
526  	
527  	         SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, NULL, NULL, &pos2, NULL);
528  	
529  	         /* only handle qterm1 == qterm here; the other case will be handled when its turn for qterm2 to be qterm */
530  	         if( qterm1 != qterm )
531  	            continue;
532  	
533  	         idx2 = quadvar2aggr[pos2];
534  	
535  	         /* variables have to belong to the same e.c. aggregation; bilinear term has to be concave */
536  	         if( idx1 >= 0 && idx2 >= 0 && idx1 == idx2 )
537  	            ++aggrnterms[idx1];
538  	         else
539  	            ++nremterms;
540  	      }
541  	   }
542  	   assert((*nlrowaggr)->nquadvars == nquadvars);
543  	
544  	   /* create all edge-concave aggregations (empty) and remaining terms */
545  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlrowaggr)->ecaggr, nfound) );
546  	   if( nremterms > 0 )
547  	   {
548  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlrowaggr)->remtermcoefs, nremterms) );
549  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlrowaggr)->remtermvars1, nremterms) );
550  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlrowaggr)->remtermvars2, nremterms) );
551  	      (*nlrowaggr)->remtermsize = nremterms;
552  	   }
553  	   (*nlrowaggr)->necaggr = nfound;
554  	
555  	   for( i = 0; i < nfound; ++i )
556  	   {
557  	      SCIP_CALL( ecaggrCreateEmpty(scip, &(*nlrowaggr)->ecaggr[i], aggrnvars[i], aggrnterms[i]) );
558  	   }
559  	
560  	   /* add quadratic variables to the edge-concave aggregations */
561  	   for( i = 0; i < nquadvars; ++i )
562  	   {
563  	      int idx;
564  	
565  	      idx = quadvar2aggr[i];
566  	
567  	      if( idx >= 0)
568  	      {
569  	         SCIP_EXPR* qterm;
570  	
571  	         SCIPdebugMsg(scip, "add quadvar %d to aggr. %d\n", i, idx);
572  	
573  	         SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, NULL, NULL, NULL, NULL);
574  	         assert(SCIPisExprVar(scip, qterm));
575  	
576  	         SCIP_CALL( ecaggrAddQuadvar((*nlrowaggr)->ecaggr[idx], SCIPgetVarExprVar(qterm)) );
577  	      }
578  	   }
579  	
580  	   /* add the bilinear/square terms to the edge-concave aggregations or in the remaining part */
581  	   for( i = 0; i < nquadvars; ++i )
582  	   {
583  	      SCIP_EXPR* qterm;
584  	      SCIP_VAR* x;
585  	      SCIP_Real coef;
586  	      int idx1;
587  	      int nadjbilin;
588  	      int* adjbilin;
589  	      int j;
590  	
591  	      SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, &coef, &nadjbilin, &adjbilin, NULL);
592  	
593  	      x = SCIPgetVarExprVar(qterm);
594  	
595  	      idx1 = quadvar2aggr[i];
596  	      if( rhsaggr )
597  	         coef = -coef;
598  	
599  	      if( idx1 >= 0 && SCIPisNegative(scip, coef) )
600  	      {
601  	         SCIP_CALL( ecaggrAddBilinTerm(scip, (*nlrowaggr)->ecaggr[idx1], x, x, coef) );
602  	         SCIPdebugMsg(scip, "add term %e *%d^2 to aggr. %d\n", coef, i, idx1);
603  	      }
604  	      else
605  	      {
606  	         SCIP_CALL( nlrowaggrAddRemBilinTerm(*nlrowaggr, x, x, coef) );
607  	         SCIPdebugMsg(scip, "add term %e *%d^2 to the remaining part\n", coef, idx1);
608  	      }
609  	
610  	      for( j = 0; j < nadjbilin; ++j )
611  	      {
612  	         SCIP_EXPR* qterm1;
613  	         SCIP_EXPR* qterm2;
614  	         int pos2;
615  	         int idx2;
616  	         SCIP_VAR* y;
617  	
618  	         SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, &qterm2, &coef, &pos2, NULL);
619  	
620  	         /* only handle qterm1 == qterm here; the other case will be handled when its turn for qterm2 to be qterm */
621  	         if( qterm1 != qterm )
622  	            continue;
623  	
624  	         y = SCIPgetVarExprVar(qterm2);
625  	
626  	         idx2 = quadvar2aggr[pos2];
627  	         if( rhsaggr )
628  	            coef = -coef;
629  	
630  	         if( idx1 >= 0 && idx2 >= 0 && idx1 == idx2 )
631  	         {
632  	            SCIP_CALL( ecaggrAddBilinTerm(scip, (*nlrowaggr)->ecaggr[idx1], x, y, coef) );
633  	            SCIPdebugMsg(scip, "add term %e *%d*%d to aggr. %d\n", coef, i, pos2, idx1);
634  	         }
635  	         else
636  	         {
637  	            SCIP_CALL( nlrowaggrAddRemBilinTerm(*nlrowaggr, x, y, coef) );
638  	            SCIPdebugMsg(scip, "add term %e *%d*%d to the remaining part\n", coef, i, pos2);
639  	         }
640  	      }
641  	   }
642  	
643  	   /* free allocated memory */
644  	   SCIPfreeBufferArray(scip, &aggrnterms);
645  	   SCIPfreeBufferArray(scip, &aggrnvars);
646  	
647  	   return SCIP_OKAY;
648  	}
649  	
650  	/** frees a nonlinear row aggregation */
651  	static
652  	SCIP_RETCODE nlrowaggrFree(
653  	   SCIP*                 scip,               /**< SCIP data structure */
654  	   SCIP_NLROWAGGR**      nlrowaggr           /**< pointer to free the nonlinear row aggregation */
655  	   )
656  	{
657  	   int i;
658  	
659  	   assert(scip != NULL);
660  	   assert(nlrowaggr != NULL);
661  	   assert(*nlrowaggr != NULL);
662  	   (*nlrowaggr)->nlrow = NULL;
663  	   assert((*nlrowaggr)->quadvars != NULL);
664  	   assert((*nlrowaggr)->nquadvars > 0);
665  	   assert((*nlrowaggr)->nremterms >= 0);
666  	
667  	   /* free remaining part */
668  	   SCIPfreeBlockMemoryArrayNull(scip, &(*nlrowaggr)->remtermcoefs, (*nlrowaggr)->remtermsize);
669  	   SCIPfreeBlockMemoryArrayNull(scip, &(*nlrowaggr)->remtermvars1, (*nlrowaggr)->remtermsize);
670  	   SCIPfreeBlockMemoryArrayNull(scip, &(*nlrowaggr)->remtermvars2, (*nlrowaggr)->remtermsize);
671  	
672  	   /* free quadratic variables */
673  	   SCIPfreeBlockMemoryArray(scip, &(*nlrowaggr)->quadvars, (*nlrowaggr)->quadvarssize);
674  	   SCIPfreeBlockMemoryArray(scip, &(*nlrowaggr)->quadvar2aggr, (*nlrowaggr)->nquadvars);
675  	
676  	   /* free linear part */
677  	   if( (*nlrowaggr)->nlinvars > 0 )
678  	   {
679  	      SCIPfreeBlockMemoryArray(scip, &(*nlrowaggr)->linvars, (*nlrowaggr)->linvarssize);
680  	      SCIPfreeBlockMemoryArray(scip, &(*nlrowaggr)->lincoefs, (*nlrowaggr)->linvarssize);
681  	   }
682  	
683  	   /* free edge-concave aggregations */
684  	   for( i = 0; i < (*nlrowaggr)->necaggr; ++i )
685  	   {
686  	      SCIP_CALL( ecaggrFree(scip, &(*nlrowaggr)->ecaggr[i]) );
687  	   }
688  	   SCIPfreeBlockMemoryArray(scip, &(*nlrowaggr)->ecaggr, (*nlrowaggr)->necaggr);
689  	
690  	   /* free nlrow aggregation */
691  	   SCIPfreeBlockMemory(scip, nlrowaggr);
692  	
693  	   return SCIP_OKAY;
694  	}
695  	
696  	#ifdef SCIP_DEBUG
697  	/** prints a nonlinear row aggregation */
698  	static
699  	void nlrowaggrPrint(
700  	   SCIP*                 scip,               /**< SCIP data structure */
701  	   SCIP_NLROWAGGR*       nlrowaggr           /**< nonlinear row aggregation */
702  	   )
703  	{
704  	   int i;
705  	
706  	   SCIPdebugMsg(scip, " nlrowaggr rhs = %e\n", nlrowaggr->rhs);
707  	   SCIPdebugMsg(scip, " #remaining terms = %d\n", nlrowaggr->nremterms);
708  	
709  	   SCIPdebugMsg(scip, "remaining terms: ");
710  	   for( i = 0; i < nlrowaggr->nremterms; ++i )
711  	      SCIPdebugMsgPrint(scip, "%e %s * %s + ", nlrowaggr->remtermcoefs[i], SCIPvarGetName(nlrowaggr->remtermvars1[i]),
712  	         SCIPvarGetName(nlrowaggr->remtermvars2[i]) );
713  	   for( i = 0; i < nlrowaggr->nlinvars; ++i )
714  	      SCIPdebugMsgPrint(scip, "%e %s + ", nlrowaggr->lincoefs[i], SCIPvarGetName(nlrowaggr->linvars[i]) );
715  	   SCIPdebugMsgPrint(scip, "\n");
716  	
717  	   for( i = 0; i < nlrowaggr->necaggr; ++i )
718  	   {
719  	      SCIPdebugMsg(scip, "print e.c. aggr %d\n", i);
720  	      ecaggrPrint(scip, nlrowaggr->ecaggr[i]);
721  	   }
722  	   return;
723  	}
724  	#endif
725  	
726  	/** creates separator data */
727  	static
728  	SCIP_RETCODE sepadataCreate(
729  	   SCIP*                 scip,               /**< SCIP data structure */
730  	   SCIP_SEPADATA**       sepadata            /**< pointer to store separator data */
731  	   )
732  	{
733  	   assert(scip != NULL);
734  	   assert(sepadata != NULL);
735  	
736  	   SCIP_CALL( SCIPallocBlockMemory(scip, sepadata) );
737  	   BMSclearMemory(*sepadata);
738  	
739  	   return SCIP_OKAY;
740  	}
741  	
742  	/** frees all nonlinear row aggregations */
743  	static
744  	SCIP_RETCODE sepadataFreeNlrows(
745  	   SCIP*                 scip,               /**< SCIP data structure */
746  	   SCIP_SEPADATA*        sepadata            /**< pointer to store separator data */
747  	   )
748  	{
749  	   assert(scip != NULL);
750  	   assert(sepadata != NULL);
751  	
752  	   /* free nonlinear row aggregations */
753  	   if( sepadata->nlrowaggrs != NULL )
754  	   {
755  	      int i;
756  	
757  	      for( i = sepadata->nnlrowaggrs - 1; i >= 0; --i )
758  	      {
759  	         SCIP_CALL( nlrowaggrFree(scip, &sepadata->nlrowaggrs[i]) );
760  	      }
761  	
762  	      SCIPfreeBlockMemoryArray(scip, &sepadata->nlrowaggrs, sepadata->nlrowaggrssize);
763  	
764  	      sepadata->nlrowaggrs = NULL;
765  	      sepadata->nnlrowaggrs = 0;
766  	      sepadata->nlrowaggrssize = 0;
767  	   }
768  	
769  	   return SCIP_OKAY;
770  	}
771  	
772  	/** frees separator data */
773  	static
774  	SCIP_RETCODE sepadataFree(
775  	   SCIP*                 scip,               /**< SCIP data structure */
776  	   SCIP_SEPADATA**       sepadata            /**< pointer to store separator data */
777  	   )
778  	{
779  	   assert(scip != NULL);
780  	   assert(sepadata != NULL);
781  	   assert(*sepadata != NULL);
782  	
783  	   /* free nonlinear row aggregations */
784  	   SCIP_CALL( sepadataFreeNlrows(scip, *sepadata) );
785  	
786  	   /* free LP interface */
787  	   if( (*sepadata)->lpi != NULL )
788  	   {
789  	      SCIP_CALL( SCIPlpiFree(&((*sepadata)->lpi)) );
790  	      (*sepadata)->lpisize = 0;
791  	   }
792  	
793  	   SCIPfreeBlockMemory(scip, sepadata);
794  	
795  	   return SCIP_OKAY;
796  	}
797  	
798  	/** adds a nonlinear row aggregation to the separator data */
799  	static
800  	SCIP_RETCODE sepadataAddNlrowaggr(
801  	   SCIP*                 scip,               /**< SCIP data structure */
802  	   SCIP_SEPADATA*        sepadata,           /**< separator data */
803  	   SCIP_NLROWAGGR*       nlrowaggr           /**< non-linear row aggregation */
804  	   )
805  	{
806  	   int i;
807  	
808  	   assert(scip != NULL);
809  	   assert(sepadata != NULL);
810  	   assert(nlrowaggr != NULL);
811  	
812  	   if( sepadata->nlrowaggrssize == 0 )
813  	   {
814  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &sepadata->nlrowaggrs, 2) ); /*lint !e506*/
815  	      sepadata->nlrowaggrssize = 2;
816  	   }
817  	   else if( sepadata->nlrowaggrssize < sepadata->nnlrowaggrs + 1 )
818  	   {
819  	      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &sepadata->nlrowaggrs, sepadata->nlrowaggrssize, 2 * sepadata->nlrowaggrssize) ); /*lint !e506 !e647*/
820  	      sepadata->nlrowaggrssize *= 2;
821  	      assert(sepadata->nlrowaggrssize >= sepadata->nnlrowaggrs + 1);
822  	   }
823  	
824  	   sepadata->nlrowaggrs[ sepadata->nnlrowaggrs ] = nlrowaggr;
825  	   ++(sepadata->nnlrowaggrs);
826  	
827  	   /* update maximum e.c. aggregation size */
828  	   for( i = 0; i < nlrowaggr->necaggr; ++i )
829  	      sepadata->maxecsize = MAX(sepadata->maxecsize, nlrowaggr->ecaggr[i]->nvars);
830  	
831  	#ifdef SCIP_STATISTIC
832  	   /* update statistics */
833  	   if( nlrowaggr->rhsaggr )
834  	      ++(sepadata->nrhsnlrowaggrs);
835  	   else
836  	      ++(sepadata->nlhsnlrowaggrs);
837  	#endif
838  	
839  	   return SCIP_OKAY;
840  	}
841  	
842  	/** returns min{val-lb,ub-val} / (ub-lb) */
843  	static
844  	SCIP_Real phi(
845  	   SCIP*                 scip,               /**< SCIP data structure */
846  	   SCIP_Real             val,                /**< solution value */
847  	   SCIP_Real             lb,                 /**< lower bound */
848  	   SCIP_Real             ub                  /**< upper bound */
849  	   )
850  	{
851  	   if( SCIPisFeasEQ(scip, lb, ub) )
852  	      return 0.0;
853  	
854  	   /* adjust */
855  	   val = MAX(val, lb);
856  	   val = MIN(val, ub);
857  	
858  	   return MIN(ub - val, val - lb) / (ub - lb);
859  	}
860  	
861  	/** creates an MIP to search for cycles with an odd number of positive edges in the graph representation of a nonlinear row
862  	 *
863  	 *  The model uses directed binary arc flow variables.
864  	 *  We introduce for all quadratic elements a forward and backward edge.
865  	 *  If the term is quadratic (e.g., loop in the graph) we fix the corresponding variables to zero.
866  	 *  This leads to an easy mapping between quadratic elements and the variables of the MIP.
867  	 */
868  	static
869  	SCIP_RETCODE createMIP(
870  	   SCIP*                 scip,               /**< SCIP data structure */
871  	   SCIP*                 subscip,            /**< auxiliary SCIP to search aggregations */
872  	   SCIP_SEPADATA*        sepadata,           /**< separator data */
873  	   SCIP_NLROW*           nlrow,              /**< nonlinear row */
874  	   SCIP_Bool             rhsaggr,            /**< consider nonlinear row aggregation for g(x) <= rhs (TRUE) or
875  	                                              *   lhs <= g(x) (FALSE) */
876  	   SCIP_VAR**            forwardarcs,        /**< array to store all forward arc variables */
877  	   SCIP_VAR**            backwardarcs,       /**< array to store all backward arc variables */
878  	   SCIP_Real*            nodeweights,        /**< weights for each node of the graph */
879  	   int*                  nedges,             /**< pointer to store the number of nonexcluded edges in the graph */
880  	   int*                  narcs               /**< pointer to store the number of created arc variables (number of square and bilinear terms) */
881  	   )
882  	{
883  	   SCIP_VAR** oddcyclearcs;
884  	   SCIP_CONS** flowcons;
885  	   SCIP_CONS* cyclelengthcons;
886  	   SCIP_CONS* oddcyclecons;
887  	   char name[SCIP_MAXSTRLEN];
888  	   SCIP_EXPR* expr;
889  	   int noddcyclearcs;
890  	   int nnodes;
891  	   int nquadexprs;
892  	   int nbilinexprs;
893  	   int i;
894  	   int arcidx;
895  	
896  	   assert(subscip != NULL);
897  	   assert(forwardarcs != NULL);
898  	   assert(backwardarcs != NULL);
899  	   assert(nedges != NULL);
900  	   assert(sepadata->minaggrsize <= sepadata->maxaggrsize);
901  	
902  	   expr = SCIPnlrowGetExpr(nlrow);
903  	   SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadexprs, &nbilinexprs, NULL, NULL);
904  	
905  	   nnodes = nquadexprs;
906  	   *nedges = 0;
907  	   *narcs = 0;
908  	
909  	   assert(nnodes > 0);
910  	
911  	   noddcyclearcs = 0;
912  	   SCIP_CALL( SCIPallocBufferArray(subscip, &oddcyclearcs, 2*nbilinexprs) );
913  	
914  	   /* create problem with default plug-ins */
915  	   SCIP_CALL( SCIPcreateProbBasic(subscip, "E.C. aggregation MIP") );
916  	   SCIP_CALL( SCIPsetObjsense(subscip, SCIP_OBJSENSE_MAXIMIZE) );
917  	   SCIP_CALL( SCIPincludeDefaultPlugins(subscip) );
918  	
919  	   /* create forward and backward arc variables */
920  	   for( i = 0; i < nquadexprs; ++i )
921  	   {
922  	      SCIP_EXPR* qterm;
923  	      SCIP_Real coef;
924  	      int nadjbilin;
925  	      int* adjbilin;
926  	      int j;
927  	
928  	      SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, &coef, &nadjbilin, &adjbilin, NULL);
929  	
930  	      if( !SCIPisZero(scip, coef) )
931  	      {
932  	         /* squares (loops) are fixed to zero */
933  	         SCIPdebugMsg(scip, "edge {%d,%d} = {%s,%s}  coeff=%e edgeweight=0\n", i, i,
934  	            SCIPvarGetName(SCIPgetVarExprVar(qterm)), SCIPvarGetName(SCIPgetVarExprVar(qterm)),
935  	            coef);
936  	
937  	         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "x#%d#%d", i, i);
938  	         SCIP_CALL( SCIPcreateVarBasic(subscip, &forwardarcs[*narcs], name, 0.0, 0.0, 0.01, SCIP_VARTYPE_BINARY) );
939  	         SCIP_CALL( SCIPaddVar(subscip, forwardarcs[*narcs]) );
940  	
941  	         SCIP_CALL( SCIPcreateVarBasic(subscip, &backwardarcs[*narcs], name, 0.0, 0.0, 0.01, SCIP_VARTYPE_BINARY) );
942  	         SCIP_CALL( SCIPaddVar(subscip, backwardarcs[*narcs]) );
943  	
944  	         ++*narcs;
945  	      }
946  	
947  	      for( j = 0 ; j < nadjbilin; ++j )
948  	      {
949  	         SCIP_EXPR* qterm1;
950  	         SCIP_EXPR* qterm2;
951  	         int pos2;
952  	         SCIP_Real edgeweight;
953  	         SCIP_CONS* noparallelcons;
954  	
955  	         SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, &qterm2, &coef, &pos2, NULL);
956  	
957  	         /* handle qterm == qterm2 later */
958  	         if( qterm1 != qterm )
959  	            continue;
960  	
961  	         edgeweight = nodeweights[i] + nodeweights[pos2];
962  	         SCIPdebugMsg(scip, "edge {%d,%d} = {%s,%s}  coeff=%e edgeweight=%e\n", i, pos2,
963  	            SCIPvarGetName(SCIPgetVarExprVar(qterm1)), SCIPvarGetName(SCIPgetVarExprVar(qterm2)),
964  	            coef, edgeweight);
965  	
966  	         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "x#%d#%d", i, pos2);
967  	         SCIP_CALL( SCIPcreateVarBasic(subscip, &forwardarcs[*narcs], name, 0.0, 1.0, 0.01 + edgeweight, SCIP_VARTYPE_BINARY) );
968  	         SCIP_CALL( SCIPaddVar(subscip, forwardarcs[*narcs]) );
969  	
970  	         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "x#%d#%d", i, pos2);
971  	         SCIP_CALL( SCIPcreateVarBasic(subscip, &backwardarcs[*narcs], name, 0.0, 1.0, 0.01 + edgeweight, SCIP_VARTYPE_BINARY) );
972  	         SCIP_CALL( SCIPaddVar(subscip, backwardarcs[*narcs]) );
973  	
974  	         ++(*nedges);
975  	
976  	         /* store all arcs which are important for the odd cycle property (no loops) */
977  	         if( rhsaggr && SCIPisPositive(scip, coef) )
978  	         {
979  	            assert(noddcyclearcs < 2*nbilinexprs-1);
980  	            oddcyclearcs[noddcyclearcs++] = forwardarcs[i];
981  	            oddcyclearcs[noddcyclearcs++] = backwardarcs[i];
982  	         }
983  	
984  	         if( !rhsaggr && SCIPisNegative(scip, coef) )
985  	         {
986  	            assert(noddcyclearcs < 2*nbilinexprs-1);
987  	            oddcyclearcs[noddcyclearcs++] = forwardarcs[i];
988  	            oddcyclearcs[noddcyclearcs++] = backwardarcs[i];
989  	         }
990  	
991  	         /* add constraints to ensure no parallel edges  */
992  	         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cons_noparalleledges");
993  	         SCIP_CALL( SCIPcreateConsBasicLinear(subscip, &noparallelcons, name, 0, NULL, NULL, 0.0, 1.0) );
994  	         SCIP_CALL( SCIPaddCoefLinear(subscip, noparallelcons, forwardarcs[*narcs], 1.0) );
995  	         SCIP_CALL( SCIPaddCoefLinear(subscip, noparallelcons, backwardarcs[*narcs], 1.0) );
996  	         SCIP_CALL( SCIPaddCons(subscip, noparallelcons) );
997  	         SCIP_CALL( SCIPreleaseCons(subscip, &noparallelcons) );
998  	
999  	         ++*narcs;
1000 	      }
1001 	   }
1002 	   assert(*narcs > 0);
1003 	
1004 	   /* odd cycle property constraint */
1005 	   (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cons_oddcycle");
1006 	   SCIP_CALL( SCIPcreateConsBasicXor(subscip, &oddcyclecons, name, TRUE, noddcyclearcs, oddcyclearcs) );
1007 	   SCIP_CALL( SCIPaddCons(subscip, oddcyclecons) );
1008 	   SCIP_CALL( SCIPreleaseCons(subscip, &oddcyclecons) );
1009 	   SCIPfreeBufferArray(subscip, &oddcyclearcs);
1010 	
1011 	   /* cycle length constraint */
1012 	   (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cons_cyclelength");
1013 	   SCIP_CALL( SCIPcreateConsBasicLinear(subscip, &cyclelengthcons, name, 0, NULL, NULL,
1014 	         (SCIP_Real) sepadata->minaggrsize, (SCIP_Real) sepadata->maxaggrsize) );
1015 	
1016 	   for( i = 0; i < *narcs; ++i )
1017 	   {
1018 	      SCIP_CALL( SCIPaddCoefLinear(subscip, cyclelengthcons, forwardarcs[i], 1.0) );
1019 	      SCIP_CALL( SCIPaddCoefLinear(subscip, cyclelengthcons, backwardarcs[i], 1.0) );
1020 	   }
1021 	
1022 	   SCIP_CALL( SCIPaddCons(subscip, cyclelengthcons) );
1023 	   SCIP_CALL( SCIPreleaseCons(subscip, &cyclelengthcons) );
1024 	
1025 	   /* create flow conservation constraints */
1026 	   SCIP_CALL( SCIPallocBufferArray(subscip, &flowcons, nnodes) );
1027 	
1028 	   for( i = 0; i < nnodes; ++i )
1029 	   {
1030 	      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cons_flowconservation#%d", i);
1031 	      SCIP_CALL( SCIPcreateConsBasicLinear(subscip, &flowcons[i], name, 0, NULL, NULL, 0.0, 0.0) );
1032 	   }
1033 	
1034 	   arcidx = 0;
1035 	   for( i = 0; i < nquadexprs; ++i )
1036 	   {
1037 	      SCIP_EXPR* qterm;
1038 	      SCIP_Real coef;
1039 	      int nadjbilin;
1040 	      int* adjbilin;
1041 	      int j;
1042 	
1043 	      SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, &coef, &nadjbilin, &adjbilin, NULL);
1044 	
1045 	      if( !SCIPisZero(scip, coef) )
1046 	         ++arcidx;
1047 	
1048 	      for( j = 0 ; j < nadjbilin; ++j )
1049 	      {
1050 	         SCIP_EXPR* qterm1;
1051 	         int pos2;
1052 	
1053 	         SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, NULL, NULL, &pos2, NULL);
1054 	
1055 	         /* handle qterm == qterm2 later */
1056 	         if( qterm1 != qterm )
1057 	            continue;
1058 	
1059 	         SCIP_CALL( SCIPaddCoefLinear(subscip, flowcons[i], forwardarcs[arcidx], 1.0) );
1060 	         SCIP_CALL( SCIPaddCoefLinear(subscip, flowcons[i], backwardarcs[arcidx], -1.0) );
1061 	
1062 	         SCIP_CALL( SCIPaddCoefLinear(subscip, flowcons[pos2], forwardarcs[arcidx], -1.0) );
1063 	         SCIP_CALL( SCIPaddCoefLinear(subscip, flowcons[pos2], backwardarcs[arcidx], 1.0) );
1064 	
1065 	         ++arcidx;
1066 	      }
1067 	   }
1068 	   assert(arcidx == *narcs);
1069 	
1070 	   for( i = 0; i < nnodes; ++i )
1071 	   {
1072 	      SCIP_CALL( SCIPaddCons(subscip, flowcons[i]) );
1073 	      SCIP_CALL( SCIPreleaseCons(subscip, &flowcons[i]) );
1074 	   }
1075 	
1076 	   SCIPfreeBufferArray(subscip, &flowcons);
1077 	
1078 	   return SCIP_OKAY;
1079 	}
1080 	
1081 	/** fixed all arc variables (u,v) for which u or v is already in an edge-concave aggregation */
1082 	static
1083 	SCIP_RETCODE updateMIP(
1084 	   SCIP*                 subscip,            /**< auxiliary SCIP to search aggregations */
1085 	   SCIP_NLROW*           nlrow,              /**< nonlinear row */
1086 	   SCIP_VAR**            forwardarcs,        /**< forward arc variables */
1087 	   SCIP_VAR**            backwardarcs,       /**< backward arc variables */
1088 	   int*                  quadvar2aggr,       /**< mapping of quadvars to e.c. aggr. index (< 0: in no aggr.) */
1089 	   int*                  nedges              /**< pointer to store the number of nonexcluded edges */
1090 	   )
1091 	{
1092 	   SCIP_EXPR* expr;
1093 	   int nquadexprs;
1094 	   int arcidx;
1095 	   int i;
1096 	
1097 	   assert(subscip != NULL);
1098 	   assert(nlrow != NULL);
1099 	   assert(forwardarcs != NULL);
1100 	   assert(backwardarcs != NULL);
1101 	   assert(quadvar2aggr != NULL);
1102 	   assert(nedges != NULL);
1103 	
1104 	   SCIP_CALL( SCIPfreeTransform(subscip) );
1105 	
1106 	   /* recompute the number of edges */
1107 	   *nedges = 0;
1108 	
1109 	   expr = SCIPnlrowGetExpr(nlrow);
1110 	   SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
1111 	
1112 	   /* fix each arc to 0 if at least one of its nodes is contained in an e.c. aggregation */
1113 	   arcidx = 0;
1114 	   for( i = 0; i < nquadexprs; ++i )
1115 	   {
1116 	      SCIP_EXPR* qterm;
1117 	      SCIP_Real coef;
1118 	      int nadjbilin;
1119 	      int* adjbilin;
1120 	      int j;
1121 	
1122 	      SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, &coef, &nadjbilin, &adjbilin, NULL);
1123 	
1124 	      if( !SCIPisZero(subscip, coef) )
1125 	      {
1126 	         if( quadvar2aggr[i] != -1 )
1127 	         {
1128 	            SCIP_CALL( SCIPchgVarUb(subscip, forwardarcs[arcidx], 0.0) );
1129 	            SCIP_CALL( SCIPchgVarUb(subscip, backwardarcs[arcidx], 0.0) );
1130 	         }
1131 	         ++arcidx;
1132 	      }
1133 	
1134 	      for( j = 0 ; j < nadjbilin; ++j )
1135 	      {
1136 	         SCIP_EXPR* qterm1;
1137 	         int pos2;
1138 	
1139 	         SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, NULL, NULL, &pos2, NULL);
1140 	
1141 	         /* handle qterm == qterm2 later */
1142 	         if( qterm1 != qterm )
1143 	            continue;
1144 	
1145 	         if( quadvar2aggr[i] != -1 || quadvar2aggr[pos2] != -1 )
1146 	         {
1147 	            SCIP_CALL( SCIPchgVarUb(subscip, forwardarcs[arcidx], 0.0) );
1148 	            SCIP_CALL( SCIPchgVarUb(subscip, backwardarcs[arcidx], 0.0) );
1149 	         }
1150 	         else
1151 	            ++*nedges;
1152 	
1153 	         ++arcidx;
1154 	      }
1155 	   }
1156 	
1157 	   return SCIP_OKAY;
1158 	}
1159 	
1160 	/** stores the best edge-concave aggregation found by the MIP model */
1161 	static
1162 	SCIP_RETCODE storeAggrFromMIP(
1163 	   SCIP*                 subscip,            /**< auxiliary SCIP to search aggregations */
1164 	   SCIP_NLROW*           nlrow,              /**< nonlinear row */
1165 	   SCIP_VAR**            forwardarcs,        /**< forward arc variables */
1166 	   SCIP_VAR**            backwardarcs,       /**< backward arc variables */
1167 	   int*                  quadvar2aggr,       /**< mapping of quadvars to e.c. aggr. index (< 0: in no aggr.) */
1168 	   int                   nfoundsofar         /**< number of e.c. aggregation found so far */
1169 	   )
1170 	{
1171 	   SCIP_SOL* sol;
1172 	   SCIP_EXPR* expr;
1173 	   int nquadexprs;
1174 	   int arcidx;
1175 	   int i;
1176 	
1177 	   assert(subscip != NULL);
1178 	   assert(nlrow != NULL);
1179 	   assert(forwardarcs != NULL);
1180 	   assert(backwardarcs != NULL);
1181 	   assert(quadvar2aggr != NULL);
1182 	   assert(nfoundsofar >= 0);
1183 	   assert(SCIPgetStatus(subscip) < SCIP_STATUS_INFEASIBLE);
1184 	   assert(SCIPgetNSols(subscip) > 0);
1185 	
1186 	   sol = SCIPgetBestSol(subscip);
1187 	   assert(sol != NULL);
1188 	
1189 	   expr = SCIPnlrowGetExpr(nlrow);
1190 	   SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
1191 	
1192 	   /* fix each arc to 0 if at least one of its nodes is contained in an e.c. aggregation */
1193 	   arcidx = 0;
1194 	   for( i = 0; i < nquadexprs; ++i )
1195 	   {
1196 	      SCIP_EXPR* qterm;
1197 	      SCIP_Real coef;
1198 	      int nadjbilin;
1199 	      int* adjbilin;
1200 	      int j;
1201 	
1202 	      SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, &coef, &nadjbilin, &adjbilin, NULL);
1203 	
1204 	      if( !SCIPisZero(subscip, coef) )
1205 	      {
1206 	         if( SCIPisGT(subscip, SCIPgetSolVal(subscip, sol, forwardarcs[arcidx]), 0.5) ||
1207 	             SCIPisGT(subscip, SCIPgetSolVal(subscip, sol, backwardarcs[arcidx]), 0.5) )
1208 	         {
1209 	            assert(quadvar2aggr[i] == -1 || quadvar2aggr[i] == nfoundsofar);
1210 	            quadvar2aggr[i] = nfoundsofar;
1211 	         }
1212 	
1213 	         ++arcidx;
1214 	      }
1215 	
1216 	      for( j = 0; j < nadjbilin; ++j )
1217 	      {
1218 	         SCIP_EXPR* qterm1;
1219 	         int pos2;
1220 	
1221 	         SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, NULL, NULL, &pos2, NULL);
1222 	
1223 	         /* handle qterm == qterm2 later */
1224 	         if( qterm1 != qterm )
1225 	            continue;
1226 	
1227 	         if( SCIPisGT(subscip, SCIPgetSolVal(subscip, sol, forwardarcs[arcidx]), 0.5) ||
1228 	             SCIPisGT(subscip, SCIPgetSolVal(subscip, sol, backwardarcs[arcidx]), 0.5) )
1229 	         {
1230 	            assert(quadvar2aggr[i] == -1 || quadvar2aggr[i] == nfoundsofar);
1231 	            assert(quadvar2aggr[pos2] == -1 || quadvar2aggr[pos2] == nfoundsofar);
1232 	
1233 	            quadvar2aggr[i] = nfoundsofar;
1234 	            quadvar2aggr[pos2] = nfoundsofar;
1235 	         }
1236 	
1237 	         ++arcidx;
1238 	      }
1239 	   }
1240 	
1241 	   return SCIP_OKAY;
1242 	}
1243 	
1244 	/** searches for edge-concave aggregations with a MIP model based on binary flow variables */
1245 	static
1246 	SCIP_RETCODE searchEcAggrWithMIP(
1247 	   SCIP*                 subscip,            /**< SCIP data structure */
1248 	   SCIP_Real             timelimit,          /**< time limit to solve the MIP */
1249 	   int                   nedges,             /**< number of nonexcluded undirected edges */
1250 	   SCIP_Bool*            aggrleft,           /**< pointer to store if there might be a left aggregation */
1251 	   SCIP_Bool*            found               /**< pointer to store if we have found an aggregation */
1252 	   )
1253 	{
1254 	   assert(subscip != NULL);
1255 	   assert(aggrleft != NULL);
1256 	   assert(found != NULL);
1257 	   assert(nedges >= 0);
1258 	
1259 	   *aggrleft = TRUE;
1260 	   *found = FALSE;
1261 	
1262 	   if( SCIPisLE(subscip, timelimit, 0.0) )
1263 	      return SCIP_OKAY;
1264 	
1265 	   /* set working limits */
1266 	   SCIP_CALL( SCIPsetRealParam(subscip, "limits/time", timelimit) );
1267 	   SCIP_CALL( SCIPsetLongintParam(subscip, "limits/totalnodes", SUBSCIP_NODELIMIT) );
1268 	
1269 	   /* set heuristics to aggressive */
1270 	   SCIP_CALL( SCIPsetHeuristics(subscip, SCIP_PARAMSETTING_AGGRESSIVE, TRUE) );
1271 	
1272 	   /* disable output to console in optimized mode, enable in SCIP's debug mode */
1273 	#ifdef SCIP_DEBUG
1274 	   SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 5) );
1275 	   SCIP_CALL( SCIPsetIntParam(subscip, "display/freq", 1) );
1276 	#else
1277 	   SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 0) );
1278 	#endif
1279 	
1280 	   SCIP_CALL( SCIPsolve(subscip) );
1281 	
1282 	   /* no more aggregation left if the MIP is infeasible */
1283 	   if( SCIPgetStatus(subscip) >= SCIP_STATUS_INFEASIBLE )
1284 	   {
1285 	      *found = FALSE;
1286 	      *aggrleft = FALSE;
1287 	      return SCIP_OKAY;
1288 	   }
1289 	
1290 	   if( SCIPgetNSols(subscip) > 0 )
1291 	   {
1292 	      *found = TRUE;
1293 	      *aggrleft = TRUE;
1294 	
1295 	#ifdef SCIP_DEBUG
1296 	      if( SCIPgetNSols(subscip) > 0 )
1297 	      {
1298 	         SCIP_CALL( SCIPprintSol(subscip, SCIPgetBestSol(subscip), NULL , FALSE) );
1299 	      }
1300 	#endif
1301 	   }
1302 	
1303 	   return SCIP_OKAY;
1304 	}
1305 	
1306 	/** creates a tclique graph from a given nonlinear row
1307 	 *
1308 	 *  SCIP's clique code can only handle integer node weights; all node weights are scaled by a factor of 100; since the
1309 	 *  clique code ignores nodes with weight of zero, we add an offset of 100 to each weight
1310 	 */
1311 	static
1312 	SCIP_RETCODE createTcliqueGraph(
1313 	   SCIP_NLROW*           nlrow,              /**< nonlinear row  */
1314 	   TCLIQUE_GRAPH**       graph,              /**< TCLIQUE graph structure */
1315 	   SCIP_Real*            nodeweights         /**< weights for each quadratic variable (nodes in the graph) */
1316 	   )
1317 	{
1318 	   SCIP_EXPR* expr;
1319 	   int nquadexprs;
1320 	   int i;
1321 	
1322 	   assert(graph != NULL);
1323 	   assert(nlrow != NULL);
1324 	
1325 	   /* create the tclique graph */
1326 	   if( !tcliqueCreate(graph) )
1327 	   {
1328 	      SCIPerrorMessage("could not create clique graph\n");
1329 	      return SCIP_ERROR;
1330 	   }
1331 	
1332 	   expr = SCIPnlrowGetExpr(nlrow);
1333 	   SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
1334 	
1335 	   /* add all nodes to the tclique graph */
1336 	   for( i = 0; i < nquadexprs; ++i )
1337 	   {
1338 	      int nodeweight;
1339 	
1340 	      /* note: clique code can only handle integer weights */
1341 	      nodeweight = 100 + (int)(100 * nodeweights[i]);
1342 	      /* SCIPdebugMsg(scip, "%d (%s): nodeweight %d \n", i, SCIPvarGetName(SCIPnlrowGetQuadVars(nlrow)[i]), nodeweight); */
1343 	
1344 	      if( !tcliqueAddNode(*graph, i, nodeweight) )
1345 	      {
1346 	         SCIPerrorMessage("could not add node to clique graph\n");
1347 	         return SCIP_ERROR;
1348 	      }
1349 	   }
1350 	
1351 	   /* add all edges */
1352 	   for( i = 0; i < nquadexprs; ++i )
1353 	   {
1354 	      SCIP_EXPR* qterm;
1355 	      int nadjbilin;
1356 	      int* adjbilin;
1357 	      int j;
1358 	
1359 	      SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, NULL, &nadjbilin, &adjbilin, NULL);
1360 	
1361 	      for( j = 0; j < nadjbilin; ++j )
1362 	      {
1363 	         SCIP_EXPR* qterm1;
1364 	         SCIP_EXPR* qterm2;
1365 	         int pos2;
1366 	
1367 	         SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, &qterm2, NULL, &pos2, NULL);
1368 	
1369 	         /* handle qterm == qterm2 later */
1370 	         if( qterm1 != qterm )
1371 	            continue;
1372 	
1373 	#ifdef SCIP_DEBUG_DETAILED
1374 	         SCIPdebugMessage("   add edge (%d, %d) = (%s,%s) to tclique graph\n",
1375 	            SCIPvarGetIndex(SCIPgetVarExprVar(qterm1)), SCIPvarGetIndex(SCIPgetVarExprVar(qterm2)),
1376 	            SCIPvarGetName(SCIPgetVarExprVar(qterm1)), SCIPvarGetName(SCIPgetVarExprVar(qterm2)));
1377 	#endif
1378 	
1379 	         if( !tcliqueAddEdge(*graph, i, pos2) )
1380 	         {
1381 	            SCIPerrorMessage("could not add edge to clique graph\n");
1382 	            return SCIP_ERROR;
1383 	         }
1384 	      }
1385 	   }
1386 	
1387 	   /* flush the clique graph */
1388 	   if( !tcliqueFlush(*graph) )
1389 	   {
1390 	      SCIPerrorMessage("could not flush the clique graph\n");
1391 	      return SCIP_ERROR;
1392 	   }
1393 	
1394 	   return SCIP_OKAY;
1395 	}
1396 	
1397 	/** searches for edge-concave aggregations by computing cliques in the graph representation of a given nonlinear row
1398 	 *
1399 	 *  update graph, compute clique, store clique; after computing a clique we heuristically check if the clique contains
1400 	 *  at least one good cycle
1401 	 */
1402 	static
1403 	SCIP_RETCODE searchEcAggrWithCliques(
1404 	   SCIP*                 scip,               /**< SCIP data structure */
1405 	   TCLIQUE_GRAPH*        graph,              /**< TCLIQUE graph structure */
1406 	   SCIP_SEPADATA*        sepadata,           /**< separator data */
1407 	   SCIP_NLROW*           nlrow,              /**< nonlinear row */
1408 	   int*                  quadvar2aggr,       /**< mapping of quadvars to e.c. aggr. index (< 0: in no aggr.) */
1409 	   int                   nfoundsofar,        /**< number of e.c. aggregation found so far */
1410 	   SCIP_Bool             rhsaggr,            /**< consider nonlinear row aggregation for g(x) <= rhs (TRUE) or
1411 	                                              *   lhs <= g(x) (FALSE) */
1412 	   SCIP_Bool*            foundaggr,          /**< pointer to store if we have found an aggregation */
1413 	   SCIP_Bool*            foundclique         /**< pointer to store if we have found a clique */
1414 	   )
1415 	{
1416 	   SCIP_HASHMAP* cliquemap;
1417 	   TCLIQUE_STATUS status;
1418 	   SCIP_EXPR* expr;
1419 	   int nquadexprs;
1420 	   int* maxcliquenodes;
1421 	   int* degrees;
1422 	   int nmaxcliquenodes;
1423 	   int maxcliqueweight;
1424 	   int noddcycleedges;
1425 	   int ntwodegrees;
1426 	   int aggrsize;
1427 	   int i;
1428 	
1429 	   assert(graph != NULL);
1430 	   assert(nfoundsofar >= 0);
1431 	   assert(foundaggr != NULL);
1432 	   assert(foundclique != NULL);
1433 	
1434 	   cliquemap = NULL;
1435 	   *foundaggr = FALSE;
1436 	   *foundclique = FALSE;
1437 	
1438 	   expr = SCIPnlrowGetExpr(nlrow);
1439 	   SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
1440 	   assert(nquadexprs == tcliqueGetNNodes(graph));
1441 	
1442 	   /* exclude all nodes which are already in an edge-concave aggregation (no flush is needed) */
1443 	   for( i = 0; i < nquadexprs; ++i )
1444 	   {
1445 	      if( quadvar2aggr[i] != -1 )
1446 	      {
1447 	         SCIPdebugMsg(scip, "exclude node %d from clique graph\n", i);
1448 	         tcliqueChangeWeight(graph, i, 0);
1449 	      }
1450 	   }
1451 	
1452 	   SCIP_CALL( SCIPallocBufferArray(scip, &maxcliquenodes, nquadexprs) );
1453 	
1454 	   /* solve clique problem */
1455 	   tcliqueMaxClique(tcliqueGetNNodes, tcliqueGetWeights, tcliqueIsEdge, tcliqueSelectAdjnodes, graph, NULL, NULL,
1456 	      maxcliquenodes, &nmaxcliquenodes, &maxcliqueweight, CLIQUE_MAXFIRSTNODEWEIGHT, CLIQUE_MINWEIGHT,
1457 	      CLIQUE_MAXNTREENODES, CLIQUE_BACKTRACKFREQ, 0, -1, NULL, &status);
1458 	
1459 	   if( status != TCLIQUE_OPTIMAL || nmaxcliquenodes < sepadata->minaggrsize )
1460 	      goto TERMINATE;
1461 	
1462 	   *foundclique = TRUE;
1463 	   aggrsize = MIN(sepadata->maxaggrsize, nmaxcliquenodes);
1464 	   SCIP_CALL( SCIPhashmapCreate(&cliquemap, SCIPblkmem(scip), aggrsize) );
1465 	
1466 	   for( i = 0; i < aggrsize; ++i )
1467 	   {
1468 	      SCIP_CALL( SCIPhashmapInsertInt(cliquemap, (void*) (size_t) maxcliquenodes[i], 0) ); /*lint !e571*/
1469 	   }
1470 	
1471 	   /* count the degree of good cycle edges for each node in the clique */
1472 	   SCIP_CALL( SCIPallocBufferArray(scip, &degrees, aggrsize) );
1473 	   BMSclearMemoryArray(degrees, aggrsize);
1474 	   ntwodegrees = 0;
1475 	
1476 	   /* count the number of positive or negative edges (depending on <= rhs or >= lhs) */
1477 	   noddcycleedges = 0;
1478 	   for( i = 0; i < nquadexprs; ++i )
1479 	   {
1480 	      SCIP_Bool isoddcycleedge;
1481 	      SCIP_EXPR* qterm;
1482 	      SCIP_Real coef;
1483 	      int nadjbilin;
1484 	      int* adjbilin;
1485 	      int j;
1486 	
1487 	      SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, &coef, &nadjbilin, &adjbilin, NULL);
1488 	
1489 	      isoddcycleedge = (rhsaggr && SCIPisPositive(scip, coef)) || (!rhsaggr && SCIPisNegative(scip, coef));
1490 	
1491 	      if( isoddcycleedge && SCIPhashmapExists(cliquemap, (void*) (size_t) i) )
1492 	      {
1493 	         ++noddcycleedges;
1494 	         ++degrees[i];
1495 	      }
1496 	
1497 	      for( j = 0; j < nadjbilin; ++j )
1498 	      {
1499 	         SCIP_EXPR* qterm1;
1500 	         SCIP_EXPR* qterm2;
1501 	         int pos2;
1502 	
1503 	         SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, &qterm2, &coef, &pos2, NULL);
1504 	
1505 	         /* handle qterm == qterm2 later */
1506 	         if( qterm1 != qterm )
1507 	            continue;
1508 	
1509 	         isoddcycleedge = (rhsaggr && SCIPisPositive(scip, coef)) || (!rhsaggr && SCIPisNegative(scip, coef));
1510 	
1511 	         if( isoddcycleedge
1512 	            && SCIPhashmapExists(cliquemap, (void*) (size_t) i)
1513 	            && SCIPhashmapExists(cliquemap, (void*) (size_t) pos2) )
1514 	         {
1515 	            ++noddcycleedges;
1516 	            ++degrees[i];
1517 	            ++degrees[pos2];
1518 	         }
1519 	      }
1520 	   }
1521 	
1522 	   /* count the number of nodes with exactly two incident odd cycle edges */
1523 	   for( i = 0; i < aggrsize; ++i )
1524 	      if( degrees[i] == 2 )
1525 	         ++ntwodegrees;
1526 	
1527 	   /* check cases for which we are sure that there are no good cycles in the clique */
1528 	   if( noddcycleedges == 0 || (aggrsize == 3 && noddcycleedges == 2) || (aggrsize == 4 && ntwodegrees == 4) )
1529 	      *foundaggr = FALSE;
1530 	   else
1531 	      *foundaggr = TRUE;
1532 	
1533 	   /* add the found clique as an edge-concave aggregation or exclude the nodes from the remaining search */
1534 	   for( i = 0; i < aggrsize; ++i )
1535 	   {
1536 	      quadvar2aggr[ maxcliquenodes[i] ] = *foundaggr ? nfoundsofar : -2;
1537 	      SCIPdebugMsg(scip, "%s %d\n", *foundaggr ? "aggregate node: " : "exclude node: ", maxcliquenodes[i]);
1538 	   }
1539 	
1540 	   SCIPfreeBufferArray(scip, &degrees);
1541 	
1542 	TERMINATE:
1543 	   if( cliquemap != NULL )
1544 	      SCIPhashmapFree(&cliquemap);
1545 	   SCIPfreeBufferArray(scip, &maxcliquenodes);
1546 	
1547 	   return SCIP_OKAY;
1548 	}
1549 	
1550 	/** helper function for searchEcAggr() */
1551 	static
1552 	SCIP_RETCODE doSeachEcAggr(
1553 	   SCIP*                 scip,               /**< SCIP data structure */
1554 	   SCIP*                 subscip,            /**< sub-SCIP data structure */
1555 	   SCIP_SEPADATA*        sepadata,           /**< separator data */
1556 	   SCIP_NLROW*           nlrow,              /**< nonlinear row */
1557 	   SCIP_SOL*             sol,                /**< current solution (might be NULL) */
1558 	   SCIP_Bool             rhsaggr,            /**< consider nonlinear row aggregation for g(x) <= rhs (TRUE) or g(x) >= lhs (FALSE) */
1559 	   int*                  quadvar2aggr,       /**< array to store for each quadratic variable in which edge-concave
1560 	                                              *   aggregation it is stored (< 0: in no aggregation); size has to be at
1561 	                                              *   least SCIPnlrowGetNQuadVars(nlrow) */
1562 	   int*                  nfound              /**< pointer to store the number of found e.c. aggregations */
1563 	   )
1564 	{
1565 	   TCLIQUE_GRAPH* graph = NULL;
1566 	   SCIP_EXPR* expr;
1567 	   SCIP_VAR** forwardarcs;
1568 	   SCIP_VAR** backwardarcs;
1569 	   SCIP_Real* nodeweights;
1570 	   SCIP_Real timelimit;
1571 	   SCIP_RETCODE retcode;
1572 	   int nunsucces = 0;
1573 	   int nedges = 0;
1574 	   int narcs;
1575 	   int nquadvars;
1576 	   int nbilinexprs;
1577 	   int i;
1578 	
1579 	   assert(subscip != NULL);
1580 	   assert(quadvar2aggr != NULL);
1581 	   assert(nfound != NULL);
1582 	
1583 	   expr = SCIPnlrowGetExpr(nlrow);
1584 	   SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadvars, &nbilinexprs, NULL, NULL);
1585 	
1586 	   retcode = SCIP_OKAY;
1587 	   *nfound = 0;
1588 	
1589 	   /* arrays to store all arc variables of the MIP model; note that we introduce variables even for loops in the graph
1590 	    * to have an easy mapping from the edges of the graph to the quadratic elements
1591 	    * nquadvars + nbilinexprs is an upper bound on the actual number of square and bilinear terms
1592 	    */
1593 	   SCIP_CALL( SCIPallocBufferArray(scip, &nodeweights, nquadvars) );
1594 	   SCIP_CALL( SCIPallocBufferArray(scip, &forwardarcs, nquadvars + nbilinexprs) );
1595 	   SCIP_CALL( SCIPallocBufferArray(scip, &backwardarcs, nquadvars + nbilinexprs) );
1596 	
1597 	   /* initialize mapping from quadvars to e.c. aggregation index (-1: quadvar is in no aggregation); compute node
1598 	    * weights
1599 	    */
1600 	   for( i = 0; i < nquadvars; ++i )
1601 	   {
1602 	      SCIP_EXPR* qterm;
1603 	      SCIP_VAR* var;
1604 	
1605 	      SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, NULL, NULL, NULL, NULL);
1606 	      assert(SCIPisExprVar(scip, qterm));
1607 	      var = SCIPgetVarExprVar(qterm);
1608 	
1609 	      quadvar2aggr[i] = -1;
1610 	      nodeweights[i] = phi(scip, SCIPgetSolVal(scip, sol, var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var));
1611 	      SCIPdebugMsg(scip, "%s = %e (%e in [%e, %e])\n", SCIPvarGetName(var), nodeweights[i], SCIPgetSolVal(scip, sol, var),
1612 	         SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var));
1613 	   }
1614 	
1615 	   SCIP_CALL( createMIP(scip, subscip, sepadata, nlrow, rhsaggr, forwardarcs, backwardarcs, nodeweights, &nedges, &narcs) );
1616 	   assert(nedges >= 0);
1617 	   assert(narcs > 0);
1618 	   SCIPdebugMsg(scip, "nedges (without loops) = %d\n", nedges);
1619 	   SCIPdebugMsg(scip, "narcs (number of quadratic terms) = %d\n", narcs);
1620 	
1621 	   SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
1622 	
1623 	   /* main loop to search for edge-concave aggregations */
1624 	   while( !SCIPisStopped(scip) )
1625 	   {
1626 	      SCIP_Bool aggrleft;
1627 	      SCIP_Bool found;
1628 	
1629 	      SCIPdebugMsg(scip, "#remaining edges = %d\n", nedges);
1630 	
1631 	      /* not enough edges left */
1632 	      if( nedges < sepadata->minaggrsize )
1633 	         break;
1634 	
1635 	      /* check whether there is enough time left; update the remaining time */
1636 	      if( !SCIPisInfinity(scip, timelimit) )
1637 	      {
1638 	         timelimit -= SCIPgetSolvingTime(scip);
1639 	         if( timelimit <= 0.0 )
1640 	         {
1641 	            SCIPdebugMsg(scip, "skip aggregation search since no time left\n");
1642 	            goto TERMINATE;
1643 	         }
1644 	      }
1645 	
1646 	      /* 1.a - search for edge-concave aggregation with the help of the MIP model */
1647 	      SCIP_CALL( searchEcAggrWithMIP(subscip, timelimit, nedges, &aggrleft, &found) );
1648 	
1649 	      /* 1.b - there are no more edge-concave aggregations left */
1650 	      if( !aggrleft )
1651 	      {
1652 	         SCIPdebugMsg(scip, "no more aggregation left\n");
1653 	         break;
1654 	      }
1655 	
1656 	      if( found )
1657 	      {
1658 	         SCIP_CALL( storeAggrFromMIP(subscip, nlrow, forwardarcs, backwardarcs, quadvar2aggr, *nfound) );
1659 	         ++(*nfound);
1660 	         nunsucces = 0;
1661 	      }
1662 	      /* try to find an edge-concave aggregation by computing cliques */
1663 	      else
1664 	      {
1665 	         SCIP_Bool foundaggr;
1666 	         SCIP_Bool foundclique;
1667 	
1668 	         ++nunsucces;
1669 	
1670 	         /* create graph if necessary */
1671 	         if( graph == NULL )
1672 	         {
1673 	            SCIP_CALL_TERMINATE( retcode, createTcliqueGraph(nlrow, &graph, nodeweights), TERMINATE );
1674 	         }
1675 	
1676 	         /* 2.a - search and store a single edge-concave aggregation by computing a clique with a good cycle */
1677 	         SCIP_CALL_FINALLY( searchEcAggrWithCliques(scip, graph, sepadata, nlrow, quadvar2aggr, *nfound, rhsaggr,
1678 	               &foundaggr, &foundclique), tcliqueFree(&graph) );
1679 	
1680 	         if( foundaggr )
1681 	         {
1682 	            assert(foundclique);
1683 	            ++(*nfound);
1684 	            nunsucces = 0;
1685 	         }
1686 	         else
1687 	            ++nunsucces;
1688 	
1689 	         /* 2.b - no clique of at least minaggrsize size found */
1690 	         if( !foundclique )
1691 	         {
1692 	            assert(!foundaggr);
1693 	            SCIPdebugMsg(scip, "did not find a clique to exclude -> leave aggregation search\n");
1694 	            break;
1695 	         }
1696 	      }
1697 	
1698 	      /* leave the algorithm if we did not find something for maxstallrounds many iterations */
1699 	      if( nunsucces >= sepadata->maxstallrounds && *nfound == 0 )
1700 	      {
1701 	         SCIPdebugMsg(scip, "did not find an e.c. aggregation for %d iterations\n", nunsucces);
1702 	         break;
1703 	      }
1704 	
1705 	      /* exclude all edges used in the last aggregation and nodes found in the clique solution */
1706 	      SCIP_CALL_FINALLY( updateMIP(subscip, nlrow, forwardarcs, backwardarcs, quadvar2aggr, &nedges), tcliqueFree(&graph) );
1707 	   }
1708 	
1709 	TERMINATE:
1710 	
1711 	#ifdef SCIP_DEBUG
1712 	   SCIPdebugMsg(scip, "aggregations found:\n");
1713 	   for( i = 0; i < nquadvars; ++i )
1714 	   {
1715 	      SCIPdebugMsg(scip, " %d in %d\n", i, quadvar2aggr[i]);
1716 	   }
1717 	#endif
1718 	
1719 	   /* free clique graph */
1720 	   if( graph != NULL )
1721 	      tcliqueFree(&graph);
1722 	
1723 	   /* free sub-SCIP */
1724 	   for( i = 0; i < narcs; ++i )
1725 	   {
1726 	      SCIP_CALL( SCIPreleaseVar(subscip, &forwardarcs[i]) );
1727 	      SCIP_CALL( SCIPreleaseVar(subscip, &backwardarcs[i]) );
1728 	   }
1729 	
1730 	   SCIPfreeBufferArray(scip, &backwardarcs);
1731 	   SCIPfreeBufferArray(scip, &forwardarcs);
1732 	   SCIPfreeBufferArray(scip, &nodeweights);
1733 	
1734 	   return retcode;
1735 	}
1736 	
1737 	/** computes a partitioning into edge-concave aggregations for a given (quadratic) nonlinear row
1738 	 *
1739 	 *  Each aggregation has to contain a cycle with an odd number of positive weighted edges (good cycles) in the corresponding graph representation.
1740 	 *  For this we use the following algorithm:
1741 	 *  -# use a MIP model based on binary flow variables to compute good cycles and store the implied subgraphs as an e.c. aggr.
1742 	 *    -# if we find a good cycle, store the implied subgraph, delete it from the graph representation and go to 1)
1743 	 *    -# if the MIP model is infeasible (there are no good cycles), STOP
1744 	 *  -# we compute a large clique C if the MIP model fails (because of working limits, etc)
1745 	 *    -# if we find a good cycle in C, store the implied subgraph of C, delete it from the graph representation and go to 1)
1746 	 *    -# if C is not large enough, STOP
1747 	 */
1748 	static
1749 	SCIP_RETCODE searchEcAggr(
1750 	   SCIP*                 scip,               /**< SCIP data structure */
1751 	   SCIP_SEPADATA*        sepadata,           /**< separator data */
1752 	   SCIP_NLROW*           nlrow,              /**< nonlinear row */
1753 	   SCIP_SOL*             sol,                /**< current solution (might be NULL) */
1754 	   SCIP_Bool             rhsaggr,            /**< consider nonlinear row aggregation for g(x) <= rhs (TRUE) or g(x) >= lhs (FALSE) */
1755 	   int*                  quadvar2aggr,       /**< array to store for each quadratic variable in which edge-concave
1756 	                                              *   aggregation it is stored (< 0: in no aggregation); size has to be at
1757 	                                              *   least SCIPnlrowGetNQuadVars(nlrow) */
1758 	   int*                  nfound              /**< pointer to store the number of found e.c. aggregations */
1759 	   )
1760 	{
1761 	   SCIP* subscip;
1762 	   SCIP_RETCODE retcode;
1763 	
1764 	   /* create and set up a sub-SCIP */
1765 	   SCIP_CALL_FINALLY( SCIPcreate(&subscip), (void)SCIPfree(&subscip) );
1766 	
1767 	   retcode = doSeachEcAggr(scip, subscip, sepadata, nlrow, sol, rhsaggr, quadvar2aggr, nfound);
1768 	
1769 	   SCIP_CALL( SCIPfree(&subscip) );
1770 	   SCIP_CALL( retcode );
1771 	
1772 	   return SCIP_OKAY;
1773 	}
1774 	
1775 	/** returns whether a given nonlinear row can be used to compute edge-concave aggregations for which their convex
1776 	 *  envelope could dominate the termwise bilinear relaxation
1777 	 *
1778 	 *  This is the case if there exists at least one cycle with
1779 	 *  an odd number of positive edges in the corresponding graph representation of the nonlinear row.
1780 	 */
1781 	static
1782 	SCIP_RETCODE isCandidate(
1783 	   SCIP*                 scip,               /**< SCIP data structure */
1784 	   SCIP_SEPADATA*        sepadata,           /**< separator data */
1785 	   SCIP_NLROW*           nlrow,              /**< nonlinear row representation of a nonlinear constraint */
1786 	   SCIP_Bool*            rhscandidate,       /**< pointer to store if we should compute edge-concave aggregations for
1787 	                                              *   the <= rhs case */
1788 	   SCIP_Bool*            lhscandidate        /**< pointer to store if we should compute edge-concave aggregations for
1789 	                                              *   the >= lhs case */
1790 	   )
1791 	{
1792 	   SCIP_EXPR* expr = NULL;
1793 	   SCIP_Bool takerow = FALSE;
1794 	   int nquadvars = 0;
1795 	   int* degrees;
1796 	   int ninterestingnodes;
1797 	   int nposedges;
1798 	   int nnegedges;
1799 	   int i;
1800 	
1801 	   assert(rhscandidate != NULL);
1802 	   assert(lhscandidate != NULL);
1803 	
1804 	   *rhscandidate = TRUE;
1805 	   *lhscandidate = TRUE;
1806 	
1807 	   /* check whether nlrow is in the NLP, is quadratic in variables, and there are enough quadratic variables */
1808 	   if( SCIPnlrowIsInNLP(nlrow) && SCIPnlrowGetExpr(nlrow) != NULL )
1809 	   {
1810 	      expr = SCIPnlrowGetExpr(nlrow);
1811 	      SCIP_CALL( SCIPcheckExprQuadratic(scip, expr, &takerow) );
1812 	   }
1813 	   if( takerow )
1814 	      takerow = SCIPexprAreQuadraticExprsVariables(expr);
1815 	   if( takerow )
1816 	   {
1817 	      SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadvars, NULL, NULL, NULL);
1818 	      takerow = nquadvars >= sepadata->minaggrsize;
1819 	   }
1820 	   if( !takerow )
1821 	   {
1822 	      *rhscandidate = FALSE;
1823 	      *lhscandidate = FALSE;
1824 	      return SCIP_OKAY;
1825 	   }
1826 	
1827 	   /* check for infinite rhs or lhs */
1828 	   if( SCIPisInfinity(scip, REALABS(SCIPnlrowGetRhs(nlrow))) )
1829 	      *rhscandidate = FALSE;
1830 	   if( SCIPisInfinity(scip, REALABS(SCIPnlrowGetLhs(nlrow))) )
1831 	      *lhscandidate = FALSE;
1832 	
1833 	   SCIP_CALL( SCIPallocClearBufferArray(scip, &degrees, nquadvars) );
1834 	
1835 	   ninterestingnodes = 0;
1836 	   nposedges = 0;
1837 	   nnegedges = 0;
1838 	
1839 	   for( i = 0; i < nquadvars; ++i )
1840 	   {
1841 	      SCIP_EXPR* qterm;
1842 	      SCIP_VAR* var1;
1843 	      int nadjbilin;
1844 	      int* adjbilin;
1845 	      int j;
1846 	
1847 	      SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, NULL, &nadjbilin, &adjbilin, NULL);
1848 	      assert(SCIPisExprVar(scip, qterm));
1849 	
1850 	      var1 = SCIPgetVarExprVar(qterm);
1851 	
1852 	      /* do not consider global fixed variables */
1853 	      if( SCIPisEQ(scip, SCIPvarGetLbGlobal(var1), SCIPvarGetUbGlobal(var1)) )
1854 	         continue;
1855 	
1856 	      for( j = 0; j < nadjbilin; ++j )
1857 	      {
1858 	         SCIP_EXPR* qterm1;
1859 	         SCIP_EXPR* qterm2;
1860 	         SCIP_VAR* var2;
1861 	         SCIP_Real coef;
1862 	         int pos2;
1863 	
1864 	         SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, &qterm2, &coef, &pos2, NULL);
1865 	
1866 	         if( qterm1 != qterm )
1867 	            continue;
1868 	
1869 	         var2 = SCIPgetVarExprVar(qterm2);
1870 	
1871 	         /* do not consider loops or global fixed variables */
1872 	         if( SCIPisEQ(scip, SCIPvarGetLbGlobal(var2), SCIPvarGetUbGlobal(var2)) )
1873 	            continue;
1874 	
1875 	         ++degrees[i];
1876 	         ++degrees[pos2];
1877 	
1878 	         /* count the number of nodes with a degree of at least 2 */
1879 	         if( degrees[i] == 2 )
1880 	            ++ninterestingnodes;
1881 	         if( degrees[pos2] == 2 )
1882 	            ++ninterestingnodes;
1883 	
1884 	         nposedges += SCIPisPositive(scip, coef) ? 1 : 0;
1885 	         nnegedges += SCIPisNegative(scip, coef) ? 1 : 0;
1886 	      }
1887 	   }
1888 	
1889 	   SCIPfreeBufferArray(scip, &degrees);
1890 	
1891 	   SCIPdebugMsg(scip, "nlrow contains: %d edges\n", nposedges + nnegedges);
1892 	
1893 	   /* too many edges, too few edges, or to few nodes with degree at least 2 in the graph  */
1894 	   if( nposedges + nnegedges > sepadata->maxbilinterms || nposedges + nnegedges < sepadata->minaggrsize
1895 	      || ninterestingnodes < sepadata->minaggrsize )
1896 	   {
1897 	      *rhscandidate = FALSE;
1898 	      *lhscandidate = FALSE;
1899 	      return SCIP_OKAY;
1900 	   }
1901 	
1902 	   /* check if there are enough positive/negative edges; for a 3-clique there has to be an odd number of those edges */
1903 	   if( nposedges == 0 || (nposedges + nnegedges == 3 && (nposedges % 2) == 0) )
1904 	      *rhscandidate = FALSE;
1905 	   if( nnegedges == 0 || (nposedges + nnegedges == 3 && (nnegedges % 2) == 0) )
1906 	      *lhscandidate = FALSE;
1907 	
1908 	   return SCIP_OKAY;
1909 	}
1910 	
1911 	/** finds and stores edge-concave aggregations for a given nonlinear row */
1912 	static
1913 	SCIP_RETCODE findAndStoreEcAggregations(
1914 	   SCIP*                 scip,               /**< SCIP data structure */
1915 	   SCIP_SEPADATA*        sepadata,           /**< separator data */
1916 	   SCIP_NLROW*           nlrow,              /**< nonlinear row */
1917 	   SCIP_SOL*             sol                 /**< current solution (might be NULL) */
1918 	   )
1919 	{
1920 	   int nquadvars;
1921 	   int* quadvar2aggr;
1922 	   SCIP_Bool rhscandidate;
1923 	   SCIP_Bool lhscandidate;
1924 	
1925 	   assert(scip != NULL);
1926 	   assert(nlrow != NULL);
1927 	   assert(sepadata != NULL);
1928 	
1929 	#ifdef SCIP_DEBUG
1930 	   SCIPdebugMsg(scip, "search for edge-concave aggregation for the nonlinear row: \n");
1931 	   SCIP_CALL( SCIPprintNlRow(scip, nlrow, NULL) );
1932 	#endif
1933 	
1934 	   /* check obvious conditions for existing cycles with an odd number of positive/negative edges */
1935 	   SCIP_CALL( isCandidate(scip, sepadata, nlrow, &rhscandidate, &lhscandidate) );
1936 	   SCIPdebugMsg(scip, "rhs candidate = %u lhs candidate = %u\n", rhscandidate, lhscandidate);
1937 	
1938 	   if( !rhscandidate && !lhscandidate )
1939 	      return SCIP_OKAY;
1940 	
1941 	   SCIPexprGetQuadraticData(SCIPnlrowGetExpr(nlrow), NULL, NULL, NULL, NULL, &nquadvars, NULL, NULL, NULL);
1942 	   SCIP_CALL( SCIPallocBufferArray(scip, &quadvar2aggr, nquadvars) ); /*lint !e705*/
1943 	
1944 	   /* search for edge-concave aggregations (consider <= rhs) */
1945 	   if( rhscandidate )
1946 	   {
1947 	      SCIP_NLROWAGGR* nlrowaggr;
1948 	      int nfound;
1949 	
1950 	      assert(!SCIPisInfinity(scip, REALABS(SCIPnlrowGetRhs(nlrow))));
1951 	
1952 	      SCIPdebugMsg(scip, "consider <= rhs\n");
1953 	      SCIP_CALL( searchEcAggr(scip, sepadata, nlrow, sol, TRUE, quadvar2aggr, &nfound) );
1954 	
1955 	      if( nfound > 0 )
1956 	      {
1957 	         SCIP_CALL( nlrowaggrCreate(scip, nlrow, &nlrowaggr, quadvar2aggr, nfound, TRUE) );
1958 	         assert(nlrow != NULL);
1959 	         SCIPdebug(nlrowaggrPrint(scip, nlrowaggr));
1960 	         SCIP_CALL( sepadataAddNlrowaggr(scip, sepadata, nlrowaggr) );
1961 	      }
1962 	   }
1963 	
1964 	   /* search for edge-concave aggregations (consider <= lhs) */
1965 	   if( lhscandidate )
1966 	   {
1967 	      SCIP_NLROWAGGR* nlrowaggr;
1968 	      int nfound;
1969 	
1970 	      assert(!SCIPisInfinity(scip, REALABS(SCIPnlrowGetLhs(nlrow))));
1971 	
1972 	      SCIPdebugMsg(scip, "consider >= lhs\n");
1973 	      SCIP_CALL( searchEcAggr(scip, sepadata, nlrow, sol, FALSE, quadvar2aggr, &nfound) );
1974 	
1975 	      if( nfound > 0 )
1976 	      {
1977 	         SCIP_CALL( nlrowaggrCreate(scip, nlrow, &nlrowaggr, quadvar2aggr, nfound, FALSE) );
1978 	         assert(nlrow != NULL);
1979 	         SCIPdebug(nlrowaggrPrint(scip, nlrowaggr));
1980 	         SCIP_CALL( sepadataAddNlrowaggr(scip, sepadata, nlrowaggr) );
1981 	      }
1982 	   }
1983 	
1984 	   SCIPfreeBufferArray(scip, &quadvar2aggr);
1985 	   return SCIP_OKAY;
1986 	}
1987 	
1988 	/*
1989 	 *  methods to compute edge-concave cuts
1990 	 */
1991 	
1992 	#ifdef SCIP_DEBUG
1993 	/** prints a given facet (candidate) */
1994 	static
1995 	void printFacet(
1996 	   SCIP*                 scip,               /**< SCIP data structure */
1997 	   SCIP_VAR**            vars,               /**< variables contained in the edge-concave aggregation */
1998 	   int                   nvars,              /**< number of variables contained in the edge-concave aggregation */
1999 	   SCIP_Real*            facet,              /**< current facet candidate */
2000 	   SCIP_Real             facetval            /**< facet evaluated at the current solution */
2001 	   )
2002 	{
2003 	   int i;
2004 	
2005 	   SCIPdebugMsg(scip, "print facet (val=%e): ", facetval);
2006 	   for( i = 0; i < nvars; ++i )
2007 	      SCIPdebugMsgPrint(scip, "%e %s + ", facet[i], SCIPvarGetName(vars[i]));
2008 	   SCIPdebugMsgPrint(scip, "%e\n", facet[nvars]);
2009 	}
2010 	#endif
2011 	
2012 	/** checks if a facet is really an underestimate for all corners of the domain [l,u]
2013 	 *
2014 	 *  Because of numerics it can happen that a facet violates a corner of the domain.
2015 	 *  To make the facet valid we subtract the maximum violation from the constant part of the facet.
2016 	 */
2017 	static
2018 	SCIP_Bool checkRikun(
2019 	   SCIP*                 scip,               /**< SCIP data structure */
2020 	   SCIP_ECAGGR*          ecaggr,             /**< edge-concave aggregation data */
2021 	   SCIP_Real*            fvals,              /**< array containing all corner values of the aggregation */
2022 	   SCIP_Real*            facet               /**< current facet candidate (of dimension ecaggr->nvars + 1) */
2023 	   )
2024 	{
2025 	   SCIP_Real maxviolation;
2026 	   SCIP_Real val;
2027 	   unsigned int i;
2028 	   unsigned int ncorner;
2029 	   unsigned int prev;
2030 	
2031 	   assert(scip != NULL);
2032 	   assert(ecaggr != NULL);
2033 	   assert(fvals != NULL);
2034 	   assert(facet != NULL);
2035 	
2036 	   ncorner = (unsigned int) poweroftwo[ecaggr->nvars];
2037 	   maxviolation = 0.0;
2038 	
2039 	   /* check for the origin */
2040 	   val = facet[ecaggr->nvars];
2041 	   for( i = 0; i < (unsigned int) ecaggr->nvars; ++i )
2042 	      val += facet[i] * SCIPvarGetLbLocal(ecaggr->vars[i]);
2043 	
2044 	   /* update  maximum violation */
2045 	   maxviolation = MAX(val - fvals[0], maxviolation);
2046 	   assert(SCIPisFeasEQ(scip, maxviolation, 0.0));
2047 	
2048 	   prev = 0;
2049 	   for( i = 1; i < ncorner; ++i )
2050 	   {
2051 	      unsigned int gray;
2052 	      unsigned int diff;
2053 	      unsigned int pos;
2054 	
2055 	      gray = i ^ (i >> 1);
2056 	      diff = gray ^ prev;
2057 	
2058 	      /* compute position of unique 1 of diff */
2059 	      pos = 0;
2060 	      while( (diff >>= 1) != 0 )
2061 	         ++pos;
2062 	
2063 	      if( gray > prev )
2064 	         val += facet[pos] * (SCIPvarGetUbLocal(ecaggr->vars[pos]) - SCIPvarGetLbLocal(ecaggr->vars[pos]));
2065 	      else
2066 	         val -= facet[pos] * (SCIPvarGetUbLocal(ecaggr->vars[pos]) - SCIPvarGetLbLocal(ecaggr->vars[pos]));
2067 	
2068 	      /* update  maximum violation */
2069 	      maxviolation = MAX(val - fvals[gray], maxviolation);
2070 	      assert(SCIPisFeasEQ(scip, maxviolation, 0.0));
2071 	
2072 	      prev = gray;
2073 	   }
2074 	
2075 	   SCIPdebugMsg(scip, "maximum violation of facet: %2.8e\n", maxviolation);
2076 	
2077 	   /* there seem to be numerical problems if the violation is too large; in this case we reject the facet */
2078 	   if( maxviolation > ADJUSTFACETTOL )
2079 	      return FALSE;
2080 	
2081 	   /* adjust constant part of the facet */
2082 	   facet[ecaggr->nvars] -= maxviolation;
2083 	
2084 	   return TRUE;
2085 	}
2086 	
2087 	/** set up LP interface to solve LPs to compute the facet of the convex envelope */
2088 	static
2089 	SCIP_RETCODE createLP(
2090 	   SCIP*                 scip,               /**< SCIP data structure */
2091 	   SCIP_SEPADATA*        sepadata            /**< separation data */
2092 	   )
2093 	{
2094 	   SCIP_Real* obj;
2095 	   SCIP_Real* lb;
2096 	   SCIP_Real* ub;
2097 	   SCIP_Real* val;
2098 	   int* beg;
2099 	   int* ind;
2100 	   int nnonz;
2101 	   int ncols;
2102 	   int nrows;
2103 	   int i;
2104 	   int k;
2105 	
2106 	   assert(scip != NULL);
2107 	   assert(sepadata != NULL);
2108 	   assert(sepadata->nnlrowaggrs > 0);
2109 	
2110 	   /* LP interface has been already created with enough rows/columns*/
2111 	   if( sepadata->lpi != NULL && sepadata->lpisize >= sepadata->maxecsize )
2112 	      return SCIP_OKAY;
2113 	
2114 	   /* size of lpi is too small; reconstruct lpi */
2115 	   if( sepadata->lpi != NULL )
2116 	   {
2117 	      SCIP_CALL( SCIPlpiFree(&sepadata->lpi) );
2118 	      sepadata->lpi = NULL;
2119 	   }
2120 	
2121 	   assert(sepadata->lpi == NULL);
2122 	   SCIP_CALL( SCIPlpiCreate(&(sepadata->lpi), SCIPgetMessagehdlr(scip), "e.c. LP", SCIP_OBJSEN_MINIMIZE) );
2123 	   sepadata->lpisize = sepadata->maxecsize;
2124 	
2125 	   nrows = sepadata->maxecsize + 1;
2126 	   ncols = poweroftwo[nrows - 1];
2127 	   nnonz = (ncols * (nrows + 1)) / 2;
2128 	   k = 0;
2129 	
2130 	   /* allocate necessary memory */
2131 	   SCIP_CALL( SCIPallocBufferArray(scip, &obj, ncols) );
2132 	   SCIP_CALL( SCIPallocBufferArray(scip, &lb, ncols) );
2133 	   SCIP_CALL( SCIPallocBufferArray(scip, &ub, ncols) );
2134 	   SCIP_CALL( SCIPallocBufferArray(scip, &beg, ncols) );
2135 	   SCIP_CALL( SCIPallocBufferArray(scip, &val, nnonz) );
2136 	   SCIP_CALL( SCIPallocBufferArray(scip, &ind, nnonz) );
2137 	
2138 	   /* calculate nonzero entries in the LP; set obj, lb, and ub to zero */
2139 	   for( i = 0; i < ncols; ++i )
2140 	   {
2141 	      int row;
2142 	      int a;
2143 	
2144 	      obj[i] = 0.0;
2145 	      lb[i] = 0.0;
2146 	      ub[i] = 0.0;
2147 	
2148 	      SCIPdebugMsg(scip, "col %i starts at position %d\n", i, k);
2149 	      beg[i] = k;
2150 	      row = 0;
2151 	      a = 1;
2152 	
2153 	      /* iterate through the bit representation of i */
2154 	      while( a <= i )
2155 	      {
2156 	         if( (a & i) != 0 )
2157 	         {
2158 	            val[k] = 1.0;
2159 	            ind[k] = row;
2160 	
2161 	            SCIPdebugMsg(scip, " val[%d][%d] = 1 (position  %d)\n", row, i, k);
2162 	
2163 	            ++k;
2164 	         }
2165 	
2166 	         a <<= 1; /*lint !e701*/
2167 	         ++row;
2168 	         assert(poweroftwo[row] == a);
2169 	      }
2170 	
2171 	      /* put 1 as a coefficient for sum_{i} \lambda_i = 1 row (last row) */
2172 	      val[k] = 1.0;
2173 	      ind[k] = nrows - 1;
2174 	      ++k;
2175 	      SCIPdebugMsg(scip, " val[%d][%d] = 1 (position  %d)\n", nrows - 1, i, k);
2176 	   }
2177 	   assert(k == nnonz);
2178 	
2179 	   /*
2180 	    * add all columns to the LP interface
2181 	    * CPLEX needs the row to exist before adding columns, so we create the rows with dummy sides
2182 	    * note that the assert is not needed once somebody fixes the LPI
2183 	    */
2184 	   assert(nrows <= ncols);
2185 	   SCIP_CALL( SCIPlpiAddRows(sepadata->lpi, nrows, obj, obj, NULL, 0, NULL, NULL, NULL) );
2186 	   SCIP_CALL( SCIPlpiAddCols(sepadata->lpi, ncols, obj, lb, ub, NULL, nnonz, beg, ind, val) );
2187 	
2188 	   /* free allocated memory */
2189 	   SCIPfreeBufferArray(scip, &ind);
2190 	   SCIPfreeBufferArray(scip, &val);
2191 	   SCIPfreeBufferArray(scip, &beg);
2192 	   SCIPfreeBufferArray(scip, &ub);
2193 	   SCIPfreeBufferArray(scip, &lb);
2194 	   SCIPfreeBufferArray(scip, &obj);
2195 	
2196 	   return SCIP_OKAY;
2197 	}
2198 	
2199 	/** evaluates an edge-concave aggregation at a corner of the domain [l,u] */
2200 	static
2201 	SCIP_Real evalCorner(
2202 	   SCIP_ECAGGR*          ecaggr,             /**< edge-concave aggregation data */
2203 	   int                   k                   /**< k-th corner */
2204 	   )
2205 	{
2206 	   SCIP_Real val;
2207 	   int i;
2208 	
2209 	   assert(ecaggr != NULL);
2210 	   assert(k >= 0 && k < poweroftwo[ecaggr->nvars]);
2211 	
2212 	   val = 0.0;
2213 	
2214 	   for( i = 0; i < ecaggr->nterms; ++i )
2215 	   {
2216 	      SCIP_Real coef;
2217 	      SCIP_Real bound1;
2218 	      SCIP_Real bound2;
2219 	      int idx1;
2220 	      int idx2;
2221 	
2222 	      idx1 = ecaggr->termvars1[i];
2223 	      idx2 = ecaggr->termvars2[i];
2224 	      coef = ecaggr->termcoefs[i];
2225 	      assert(idx1 >= 0 && idx1 < ecaggr->nvars);
2226 	      assert(idx2 >= 0 && idx2 < ecaggr->nvars);
2227 	
2228 	      bound1 = ((poweroftwo[idx1]) & k) == 0 ? SCIPvarGetLbLocal(ecaggr->vars[idx1]) : SCIPvarGetUbLocal(ecaggr->vars[idx1]); /*lint !e661*/
2229 	      bound2 = ((poweroftwo[idx2]) & k) == 0 ? SCIPvarGetLbLocal(ecaggr->vars[idx2]) : SCIPvarGetUbLocal(ecaggr->vars[idx2]); /*lint !e661*/
2230 	
2231 	      val += coef * bound1 * bound2;
2232 	   }
2233 	
2234 	   return val;
2235 	}
2236 	
2237 	/** returns (val - lb) / (ub - lb) for a in [lb, ub] */
2238 	static
2239 	SCIP_Real transformValue(
2240 	   SCIP*                 scip,               /**< SCIP data structure */
2241 	   SCIP_Real             lb,                 /**< lower bound */
2242 	   SCIP_Real             ub,                 /**< upper bound */
2243 	   SCIP_Real             val                 /**< value in [lb,ub] */
2244 	   )
2245 	{
2246 	   assert(scip != NULL);
2247 	   assert(!SCIPisInfinity(scip, -lb));
2248 	   assert(!SCIPisInfinity(scip, ub));
2249 	   assert(!SCIPisInfinity(scip, REALABS(val)));
2250 	   assert(!SCIPisFeasEQ(scip, ub - lb, 0.0)); /* this would mean that a variable has been fixed */
2251 	
2252 	   /* adjust val */
2253 	   val = MIN(val, ub);
2254 	   val = MAX(val, lb);
2255 	
2256 	   val = (val - lb) / (ub - lb);
2257 	   assert(val >= 0.0 && val <= 1.0);
2258 	
2259 	   return val;
2260 	}
2261 	
2262 	/** computes a facet of the convex envelope of an edge concave aggregation
2263 	 *
2264 	 *  The algorithm solves the following LP:
2265 	 *  \f{align}{
2266 	 *      \min  & \sum_i \lambda_i f(v_i)\\
2267 	 *      s.t.  & \sum_i \lambda_i v_i = x\\
2268 	 *            & \sum_i \lambda_i     = 1\\
2269 	 *            & \lambda \geq 0
2270 	 *  \f}
2271 	 *  where \f$f\f$ is an edge concave function, \f$x\in [l,u]\f$ is a solution of the current relaxation, and \f$v_i\f$ are the vertices of \f$[l,u]\f$.
2272 	 *  The method transforms the problem to the domain \f$[0,1]^n\f$, computes a facet, and transforms this facet to the
2273 	 *  original space. The dual solution of the LP above are the coefficients of the facet.
2274 	 *
2275 	 *  The complete algorithm works as follows:
2276 	 *  -# compute \f$f(v_i)\f$ for each corner \f$v_i\f$ of \f$[l,u]\f$
2277 	 *  -# set up the described LP for the transformed space
2278 	 *  -# solve the LP and store the resulting facet for the transformed space
2279 	 *  -# transform the facet to original space
2280 	 *  -# adjust and check facet with the algorithm of Rikun et al.
2281 	 */
2282 	static
2283 	SCIP_RETCODE computeConvexEnvelopeFacet(
2284 	   SCIP*                 scip,               /**< SCIP data structure */
2285 	   SCIP_SEPADATA*        sepadata,           /**< separation data */
2286 	   SCIP_SOL*             sol,                /**< solution (might be NULL) */
2287 	   SCIP_ECAGGR*          ecaggr,             /**< edge-concave aggregation data */
2288 	   SCIP_Real*            facet,              /**< array to store the coefficients of the resulting facet; size has to be at least (ecaggr->nvars + 1) */
2289 	   SCIP_Real*            facetval,           /**< pointer to store the value of the facet evaluated at the current solution */
2290 	   SCIP_Bool*            success             /**< pointer to store if we have found a facet */
2291 	   )
2292 	{
2293 	   SCIP_Real* fvals;
2294 	   SCIP_Real* side;
2295 	   SCIP_Real* lb;
2296 	   SCIP_Real* ub;
2297 	   SCIP_Real perturbation;
2298 	   int* inds;
2299 	   int ncorner;
2300 	   int ncols;
2301 	   int nrows;
2302 	   int i;
2303 	
2304 	   assert(scip != NULL);
2305 	   assert(sepadata != NULL);
2306 	   assert(ecaggr != NULL);
2307 	   assert(facet != NULL);
2308 	   assert(facetval != NULL);
2309 	   assert(success != NULL);
2310 	   assert(ecaggr->nvars <= sepadata->maxecsize);
2311 	
2312 	   *facetval = -SCIPinfinity(scip);
2313 	   *success = FALSE;
2314 	
2315 	   /* create LP if this has not been done yet */
2316 	   SCIP_CALL( createLP(scip, sepadata) );
2317 	
2318 	   assert(sepadata->lpi != NULL);
2319 	   assert(sepadata->lpisize >= ecaggr->nvars);
2320 	
2321 	   SCIP_CALL( SCIPlpiGetNCols(sepadata->lpi, &ncols) );
2322 	   SCIP_CALL( SCIPlpiGetNRows(sepadata->lpi, &nrows) );
2323 	   ncorner = poweroftwo[ecaggr->nvars];
2324 	
2325 	   assert(ncorner <= ncols);
2326 	   assert(ecaggr->nvars + 1 <= nrows);
2327 	   assert(nrows <= ncols);
2328 	
2329 	   /* allocate necessary memory */
2330 	   SCIP_CALL( SCIPallocBufferArray(scip, &fvals, ncols) );
2331 	   SCIP_CALL( SCIPallocBufferArray(scip, &inds, ncols) );
2332 	   SCIP_CALL( SCIPallocBufferArray(scip, &lb, ncols) );
2333 	   SCIP_CALL( SCIPallocBufferArray(scip, &ub, ncols) );
2334 	   SCIP_CALL( SCIPallocBufferArray(scip, &side, ncols) );
2335 	
2336 	   /*
2337 	    *  1. compute f(v_i) for each corner v_i of [l,u]
2338 	    *  2. set up the described LP for the transformed space
2339 	    */
2340 	   for( i = 0; i < ncols; ++i )
2341 	   {
2342 	      fvals[i] = i < ncorner ? evalCorner(ecaggr, i) : 0.0;
2343 	      inds[i] = i;
2344 	
2345 	      /* update bounds; fix variables to zero which are currently not in the LP */
2346 	      lb[i] = 0.0;
2347 	      ub[i] = i < ncorner ? 1.0 : 0.0;
2348 	      SCIPdebugMsg(scip, "bounds of LP col %d = [%e, %e]; obj = %e\n", i, lb[i], ub[i], fvals[i]);
2349 	   }
2350 	
2351 	   /* update lhs and rhs */
2352 	   perturbation = 0.001;
2353 	   for( i = 0; i < nrows; ++i )
2354 	   {
2355 	      /* note that the last row corresponds to sum_{j} \lambda_j = 1 */
2356 	      if( i < ecaggr->nvars )
2357 	      {
2358 	         SCIP_VAR* x;
2359 	
2360 	         x = ecaggr->vars[i];
2361 	         assert(x != NULL);
2362 	
2363 	         side[i] = transformValue(scip, SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x), SCIPgetSolVal(scip, sol, x));
2364 	
2365 	         /* perturb point to enforce an LP solution with ecaggr->nvars + 1 nonzero */
2366 	         side[i] += side[i] > perturbation ? -perturbation : perturbation;
2367 	         perturbation /= 1.2;
2368 	      }
2369 	      else
2370 	      {
2371 	         side[i] = (i == nrows - 1) ? 1.0 : 0.0;
2372 	      }
2373 	
2374 	      SCIPdebugMsg(scip, "LP row %d in [%e, %e]\n", i, side[i], side[i]);
2375 	   }
2376 	
2377 	   /* update LP */
2378 	   SCIP_CALL( SCIPlpiChgObj(sepadata->lpi, ncols, inds, fvals) );
2379 	   SCIP_CALL( SCIPlpiChgBounds(sepadata->lpi, ncols, inds, lb, ub) );
2380 	   SCIP_CALL( SCIPlpiChgSides(sepadata->lpi, nrows, inds, side, side) );
2381 	
2382 	   /* free memory used to build the LP */
2383 	   SCIPfreeBufferArray(scip, &side);
2384 	   SCIPfreeBufferArray(scip, &ub);
2385 	   SCIPfreeBufferArray(scip, &lb);
2386 	   SCIPfreeBufferArray(scip, &inds);
2387 	
2388 	   /*
2389 	    *  3. solve the LP and store the resulting facet for the transformed space
2390 	    */
2391 	   if( USEDUALSIMPLEX ) /*lint !e774 !e506*/
2392 	   {
2393 	      SCIP_CALL( SCIPlpiSolveDual(sepadata->lpi) );
2394 	   }
2395 	   else
2396 	   {
2397 	      SCIP_CALL( SCIPlpiSolvePrimal(sepadata->lpi) );
2398 	   }
2399 	
2400 	   /* the dual solution corresponds to the coefficients of the facet in the transformed problem; note that it might be
2401 	    * the case that the dual solution has more components than the facet array
2402 	    */
2403 	   if( ecaggr->nvars + 1 == ncols )
2404 	   {
2405 	      SCIP_CALL( SCIPlpiGetSol(sepadata->lpi, NULL, NULL, facet, NULL, NULL) );
2406 	   }
2407 	   else
2408 	   {
2409 	      SCIP_Real* dualsol;
2410 	
2411 	      SCIP_CALL( SCIPallocBufferArray(scip, &dualsol, nrows) );
2412 	
2413 	      /* get the dual solution */
2414 	      SCIP_CALL( SCIPlpiGetSol(sepadata->lpi, NULL, NULL, dualsol, NULL, NULL) );
2415 	
2416 	      for( i = 0; i < ecaggr->nvars; ++i )
2417 	         facet[i] = dualsol[i];
2418 	
2419 	      /* constant part of the facet is the last component of the dual solution */
2420 	      facet[ecaggr->nvars] = dualsol[nrows - 1];
2421 	
2422 	      SCIPfreeBufferArray(scip, &dualsol);
2423 	   }
2424 	
2425 	#ifdef SCIP_DEBUG
2426 	   SCIPdebugMsg(scip, "facet for the transformed problem: ");
2427 	   for( i = 0; i < ecaggr->nvars; ++i )
2428 	   {
2429 	      SCIPdebugMsgPrint(scip, "%3.4e * %s + ", facet[i], SCIPvarGetName(ecaggr->vars[i]));
2430 	   }
2431 	   SCIPdebugMsgPrint(scip, "%3.4e\n", facet[ecaggr->nvars]);
2432 	#endif
2433 	
2434 	   /*
2435 	    *  4. transform the facet to original space
2436 	    *  we now have the linear underestimator L(x) = beta^T x + beta_0, which needs to be transform to the original space
2437 	    *  the underestimator in the original space, G(x) = alpha^T x + alpha_0, is given by G(x) = L(T(x)), where T(.) is
2438 	    *  the transformation applied in step 2; therefore,
2439 	    *  alpha_i = beta_i/(ub_i - lb_i)
2440 	    *  alpha_0 = beta_0 - sum_i lb_i * beta_i/(ub_i - lb_i)
2441 	    */
2442 	
2443 	   SCIPdebugMsg(scip, "facet in orig. space: ");
2444 	   *facetval = 0.0;
2445 	
2446 	   for( i = 0; i < ecaggr->nvars; ++i )
2447 	   {
2448 	      SCIP_Real varlb;
2449 	      SCIP_Real varub;
2450 	
2451 	      varlb = SCIPvarGetLbLocal(ecaggr->vars[i]);
2452 	      varub = SCIPvarGetUbLocal(ecaggr->vars[i]);
2453 	      assert(!SCIPisEQ(scip, varlb, varub));
2454 	
2455 	      /* substract (\beta_i * lb_i) / (ub_i - lb_i) from current alpha_0 */
2456 	      facet[ecaggr->nvars] -= (facet[i] * varlb) / (varub - varlb);
2457 	
2458 	      /* set \alpha_i := \beta_i / (ub_i - lb_i) */
2459 	      facet[i] = facet[i] / (varub - varlb);
2460 	      *facetval += facet[i] * SCIPgetSolVal(scip, sol, ecaggr->vars[i]);
2461 	
2462 	      SCIPdebugMsgPrint(scip, "%3.4e * %s + ", facet[i], SCIPvarGetName(ecaggr->vars[i]));
2463 	   }
2464 	
2465 	   /* add constant part to the facet value */
2466 	   *facetval += facet[ecaggr->nvars];
2467 	   SCIPdebugMsgPrint(scip, "%3.4e\n", facet[ecaggr->nvars]);
2468 	
2469 	   /*
2470 	    *  5. adjust and check facet with the algorithm of Rikun et al.
2471 	    */
2472 	
2473 	   if( checkRikun(scip, ecaggr, fvals, facet) )
2474 	   {
2475 	      SCIPdebugMsg(scip, "facet pass the check of Rikun et al.\n");
2476 	      *success = TRUE;
2477 	   }
2478 	
2479 	   /* free allocated memory */
2480 	   SCIPfreeBufferArray(scip, &fvals);
2481 	
2482 	   return SCIP_OKAY;
2483 	}
2484 	
2485 	/*
2486 	 * miscellaneous methods
2487 	 */
2488 	
2489 	/** method to add a facet of the convex envelope of an edge-concave aggregation to a given cut */
2490 	static
2491 	SCIP_RETCODE addFacetToCut(
2492 	   SCIP*                 scip,               /**< SCIP data structure */
2493 	   SCIP_SOL*             sol,                /**< current solution (might be NULL) */
2494 	   SCIP_ROW*             cut,                /**< current cut (modifiable) */
2495 	   SCIP_Real*            facet,              /**< coefficient of the facet (dimension nvars + 1) */
2496 	   SCIP_VAR**            vars,               /**< variables of the facet */
2497 	   int                   nvars,              /**< number of variables in the facet */
2498 	   SCIP_Real*            cutconstant,        /**< pointer to update the constant part of the facet */
2499 	   SCIP_Real*            cutactivity,        /**< pointer to update the activity of the cut */
2500 	   SCIP_Bool*            success             /**< pointer to store if everything went fine */
2501 	   )
2502 	{
2503 	   int i;
2504 	
2505 	   assert(cut != NULL);
2506 	   assert(facet != NULL);
2507 	   assert(vars != NULL);
2508 	   assert(nvars > 0);
2509 	   assert(cutconstant != NULL);
2510 	   assert(cutactivity != NULL);
2511 	   assert(success != NULL);
2512 	
2513 	   *success = TRUE;
2514 	
2515 	   for( i = 0; i < nvars; ++i )
2516 	   {
2517 	      if( SCIPisInfinity(scip, REALABS(facet[i])) )
2518 	      {
2519 	         *success = FALSE;
2520 	         return SCIP_OKAY;
2521 	      }
2522 	
2523 	      if( !SCIPisZero(scip, facet[i]) )
2524 	      {
2525 	         /* add only a constant if the variable has been fixed */
2526 	         if( SCIPvarGetLbLocal(vars[i]) == SCIPvarGetUbLocal(vars[i]) ) /*lint !e777*/
2527 	         {
2528 	            assert(SCIPisFeasEQ(scip, SCIPvarGetLbLocal(vars[i]), SCIPgetSolVal(scip, sol, vars[i])));
2529 	            *cutconstant += facet[i] * SCIPgetSolVal(scip, sol, vars[i]);
2530 	            *cutactivity += facet[i] * SCIPgetSolVal(scip, sol, vars[i]);
2531 	         }
2532 	         else
2533 	         {
2534 	            *cutactivity += facet[i] * SCIPgetSolVal(scip, sol, vars[i]);
2535 	            SCIP_CALL( SCIPaddVarToRow(scip, cut, vars[i], facet[i]) );
2536 	         }
2537 	      }
2538 	   }
2539 	
2540 	   /* add constant part of the facet */
2541 	   *cutconstant += facet[nvars];
2542 	   *cutactivity += facet[nvars];
2543 	
2544 	   return SCIP_OKAY;
2545 	}
2546 	
2547 	/** method to add a linear term to a given cut */
2548 	static
2549 	SCIP_RETCODE addLinearTermToCut(
2550 	   SCIP*                 scip,               /**< SCIP data structure */
2551 	   SCIP_SOL*             sol,                /**< current solution (might be NULL) */
2552 	   SCIP_ROW*             cut,                /**< current cut (modifiable) */
2553 	   SCIP_VAR*             x,                  /**< linear variable */
2554 	   SCIP_Real             coeff,              /**< coefficient */
2555 	   SCIP_Real*            cutconstant,        /**< pointer to update the constant part of the facet */
2556 	   SCIP_Real*            cutactivity,        /**< pointer to update the activity of the cut */
2557 	   SCIP_Bool*            success             /**< pointer to store if everything went fine */
2558 	   )
2559 	{
2560 	   SCIP_Real activity;
2561 	
2562 	   assert(cut != NULL);
2563 	   assert(x != NULL);
2564 	   assert(!SCIPisZero(scip, coeff));
2565 	   assert(!SCIPisInfinity(scip, coeff));
2566 	   assert(cutconstant != NULL);
2567 	   assert(cutactivity != NULL);
2568 	   assert(success != NULL);
2569 	
2570 	   *success = TRUE;
2571 	   activity = SCIPgetSolVal(scip, sol, x) * coeff;
2572 	
2573 	   /* do not add a term if the activity is -infinity */
2574 	   if( SCIPisInfinity(scip, -1.0 * REALABS(activity)) )
2575 	   {
2576 	      *success = FALSE;
2577 	      return SCIP_OKAY;
2578 	   }
2579 	
2580 	   /* add activity to the constant part if the variable has been fixed */
2581 	   if( SCIPvarGetLbLocal(x) == SCIPvarGetUbLocal(x) ) /*lint !e777*/
2582 	   {
2583 	      assert(SCIPisFeasEQ(scip, SCIPvarGetLbLocal(x), SCIPgetSolVal(scip, sol, x)));
2584 	      *cutconstant += activity;
2585 	      SCIPdebugMsg(scip, "add to cut: %e\n", activity);
2586 	   }
2587 	   else
2588 	   {
2589 	      SCIP_CALL( SCIPaddVarToRow(scip, cut, x, coeff) );
2590 	      SCIPdebugMsg(scip, "add to cut: %e * %s\n", coeff, SCIPvarGetName(x));
2591 	   }
2592 	
2593 	   *cutactivity += activity;
2594 	
2595 	   return SCIP_OKAY;
2596 	}
2597 	
2598 	/** method to add an underestimate of a bilinear term to a given cut */
2599 	static
2600 	SCIP_RETCODE addBilinearTermToCut(
2601 	   SCIP*                 scip,               /**< SCIP data structure */
2602 	   SCIP_SOL*             sol,                /**< current solution (might be NULL) */
2603 	   SCIP_ROW*             cut,                /**< current cut (modifiable) */
2604 	   SCIP_VAR*             x,                  /**< first bilinear variable */
2605 	   SCIP_VAR*             y,                  /**< seconds bilinear variable */
2606 	   SCIP_Real             coeff,              /**< coefficient */
2607 	   SCIP_Real*            cutconstant,        /**< pointer to update the constant part of the facet */
2608 	   SCIP_Real*            cutactivity,        /**< pointer to update the activity of the cut */
2609 	   SCIP_Bool*            success             /**< pointer to store if everything went fine */
2610 	   )
2611 	{
2612 	   SCIP_Real activity;
2613 	
2614 	   assert(cut != NULL);
2615 	   assert(x != NULL);
2616 	   assert(y != NULL);
2617 	   assert(!SCIPisZero(scip, coeff));
2618 	   assert(cutconstant != NULL);
2619 	   assert(cutactivity != NULL);
2620 	   assert(success != NULL);
2621 	
2622 	   *success = TRUE;
2623 	   activity = coeff * SCIPgetSolVal(scip, sol, x) * SCIPgetSolVal(scip, sol, y);
2624 	
2625 	   if( SCIPisInfinity(scip, REALABS(coeff)) )
2626 	   {
2627 	      *success = FALSE;
2628 	      return SCIP_OKAY;
2629 	   }
2630 	
2631 	   /* do not add a term if the activity is -infinity */
2632 	   if( SCIPisInfinity(scip, -1.0 * REALABS(activity)) )
2633 	   {
2634 	      *success = FALSE;
2635 	      return SCIP_OKAY;
2636 	   }
2637 	
2638 	   /* quadratic case */
2639 	   if( x == y )
2640 	   {
2641 	      SCIP_Real refpoint;
2642 	      SCIP_Real lincoef;
2643 	      SCIP_Real linconst;
2644 	
2645 	      lincoef = 0.0;
2646 	      linconst = 0.0;
2647 	      refpoint = SCIPgetSolVal(scip, sol, x);
2648 	
2649 	      /* adjust the reference point */
2650 	      refpoint = SCIPisLT(scip, refpoint, SCIPvarGetLbLocal(x)) ? SCIPvarGetLbLocal(x) : refpoint;
2651 	      refpoint = SCIPisGT(scip, refpoint, SCIPvarGetUbLocal(x)) ? SCIPvarGetUbLocal(x) : refpoint;
2652 	      assert(SCIPisLE(scip, refpoint, SCIPvarGetUbLocal(x)) && SCIPisGE(scip, refpoint, SCIPvarGetLbLocal(x)));
2653 	
2654 	      if( SCIPisPositive(scip, coeff) )
2655 	         SCIPaddSquareLinearization(scip, coeff, refpoint, SCIPvarIsIntegral(x), &lincoef, &linconst, success);
2656 	      else
2657 	         SCIPaddSquareSecant(scip, coeff, SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x), &lincoef, &linconst, success);
2658 	
2659 	      *cutactivity += lincoef * refpoint + linconst;
2660 	      *cutconstant += linconst;
2661 	
2662 	      /* add underestimate to cut */
2663 	      SCIP_CALL( SCIPaddVarToRow(scip, cut, x, lincoef) );
2664 	
2665 	      SCIPdebugMsg(scip, "add to cut: %e * %s + %e\n", lincoef, SCIPvarGetName(x), linconst);
2666 	   }
2667 	   /* bilinear case */
2668 	   else
2669 	   {
2670 	      SCIP_Real refpointx;
2671 	      SCIP_Real refpointy;
2672 	      SCIP_Real lincoefx;
2673 	      SCIP_Real lincoefy;
2674 	      SCIP_Real linconst;
2675 	
2676 	      lincoefx = 0.0;
2677 	      lincoefy = 0.0;
2678 	      linconst = 0.0;
2679 	      refpointx = SCIPgetSolVal(scip, sol, x);
2680 	      refpointy = SCIPgetSolVal(scip, sol, y);
2681 	
2682 	      /* adjust the reference points */
2683 	      refpointx = SCIPisLT(scip, refpointx, SCIPvarGetLbLocal(x)) ? SCIPvarGetLbLocal(x) : refpointx;
2684 	      refpointx = SCIPisGT(scip, refpointx, SCIPvarGetUbLocal(x)) ? SCIPvarGetUbLocal(x) : refpointx;
2685 	      refpointy = SCIPisLT(scip, refpointy, SCIPvarGetLbLocal(y)) ? SCIPvarGetLbLocal(y) : refpointy;
2686 	      refpointy = SCIPisGT(scip, refpointy, SCIPvarGetUbLocal(y)) ? SCIPvarGetUbLocal(y) : refpointy;
2687 	      assert(SCIPisLE(scip, refpointx, SCIPvarGetUbLocal(x)) && SCIPisGE(scip, refpointx, SCIPvarGetLbLocal(x)));
2688 	      assert(SCIPisLE(scip, refpointy, SCIPvarGetUbLocal(y)) && SCIPisGE(scip, refpointy, SCIPvarGetLbLocal(y)));
2689 	
2690 	      SCIPaddBilinMcCormick(scip, coeff, SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x), refpointx, SCIPvarGetLbLocal(y),
2691 	         SCIPvarGetUbLocal(y), refpointy, FALSE, &lincoefx, &lincoefy, &linconst, success);
2692 	
2693 	      *cutactivity += lincoefx * refpointx + lincoefy * refpointy + linconst;
2694 	      *cutconstant += linconst;
2695 	
2696 	      /* add underestimate to cut */
2697 	      SCIP_CALL( SCIPaddVarToRow(scip, cut, x, lincoefx) );
2698 	      SCIP_CALL( SCIPaddVarToRow(scip, cut, y, lincoefy) );
2699 	
2700 	      SCIPdebugMsg(scip, "add to cut: %e * %s + %e * %s + %e\n", lincoefx, SCIPvarGetName(x), lincoefy,
2701 	         SCIPvarGetName(y), linconst);
2702 	   }
2703 	
2704 	   return SCIP_OKAY;
2705 	}
2706 	
2707 	/** method to compute and add a cut for a nonlinear row aggregation and a given solution
2708 	 *
2709 	 *  we compute for each edge concave aggregation one facet;
2710 	 *  the remaining bilinear terms will be underestimated with McCormick, secants or linearizations;
2711 	 *  constant and linear terms will be added to the cut directly
2712 	 */
2713 	static
2714 	SCIP_RETCODE computeCut(
2715 	   SCIP*                 scip,               /**< SCIP data structure */
2716 	   SCIP_SEPA*            sepa,               /**< separator */
2717 	   SCIP_SEPADATA*        sepadata,           /**< separator data */
2718 	   SCIP_NLROWAGGR*       nlrowaggr,          /**< nonlinear row aggregation */
2719 	   SCIP_SOL*             sol,                /**< current solution (might be NULL) */
2720 	   SCIP_Bool*            separated,          /**< pointer to store if we could separate the current solution */
2721 	   SCIP_Bool*            cutoff              /**< pointer to store if the current node gets cut off */
2722 	   )
2723 	{
2724 	   SCIP_ROW* cut;
2725 	   SCIP_Real* bestfacet;
2726 	   SCIP_Real bestfacetval;
2727 	   SCIP_Real cutconstant;
2728 	   SCIP_Real cutactivity;
2729 	   int bestfacetsize;
2730 	   char cutname[SCIP_MAXSTRLEN];
2731 	   SCIP_Bool found;
2732 	   SCIP_Bool islocalcut;
2733 	   int i;
2734 	
2735 	   assert(separated != NULL);
2736 	   assert(cutoff != NULL);
2737 	   assert(nlrowaggr->necaggr > 0);
2738 	   assert(nlrowaggr->nlrow != NULL);
2739 	   assert(SCIPnlrowIsInNLP(nlrowaggr->nlrow));
2740 	
2741 	   *separated = FALSE;
2742 	   *cutoff = FALSE;
2743 	   /* we use SCIPgetDepth because we add the cut to the global cut pool if cut is globally valid */
2744 	   islocalcut = SCIPgetDepth(scip) != 0;
2745 	
2746 	   /* create the cut */
2747 	   (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "ec");
2748 	   SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, sepa, cutname, -SCIPinfinity(scip), SCIPinfinity(scip), islocalcut, FALSE,
2749 	         sepadata->dynamiccuts) );
2750 	   SCIP_CALL( SCIPcacheRowExtensions(scip, cut) );
2751 	
2752 	   /* track rhs and activity of the cut */
2753 	   cutconstant = nlrowaggr->constant;
2754 	   cutactivity = 0.0;
2755 	
2756 	   /* allocate necessary memory */
2757 	   bestfacetsize = sepadata->maxaggrsize + 1;
2758 	   SCIP_CALL( SCIPallocBufferArray(scip, &bestfacet, bestfacetsize) );
2759 	
2760 	#ifdef SCIP_DEBUG
2761 	   SCIP_CALL( SCIPprintNlRow(scip, nlrowaggr->nlrow, NULL) );
2762 	
2763 	   SCIPdebugMsg(scip, "current solution:\n");
2764 	   for( i = 0; i < SCIPgetNVars(scip); ++i )
2765 	   {
2766 	      SCIP_VAR* var = SCIPgetVars(scip)[i];
2767 	      SCIPdebugMsg(scip, "  %s = [%e, %e] solval = %e\n", SCIPvarGetName(var), SCIPvarGetLbLocal(var),
2768 	         SCIPvarGetUbLocal(var), SCIPgetSolVal(scip, sol, var));
2769 	   }
2770 	#endif
2771 	
2772 	   /* compute a facet for each edge-concave aggregation */
2773 	   for( i = 0; i < nlrowaggr->necaggr; ++i )
2774 	   {
2775 	      SCIP_ECAGGR* ecaggr;
2776 	      SCIP_Bool success;
2777 	
2778 	      ecaggr = nlrowaggr->ecaggr[i];
2779 	      assert(ecaggr != NULL);
2780 	
2781 	      /* compute a facet of the convex envelope */
2782 	      SCIP_CALL( computeConvexEnvelopeFacet(scip, sepadata, sol, ecaggr, bestfacet, &bestfacetval, &found) );
2783 	
2784 	      SCIPdebugMsg(scip, "found facet for edge-concave aggregation %d/%d ? %s\n", i, nlrowaggr->necaggr,
2785 	         found ? "yes" : "no");
2786 	
2787 	#ifdef SCIP_DEBUG
2788 	      if( found )
2789 	         printFacet(scip, ecaggr->vars, ecaggr->nvars, bestfacet, bestfacetval);
2790 	#endif
2791 	
2792 	      /* do not add any cut because we did not found a facet for at least one edge-concave aggregation */
2793 	      if( !found ) /*lint !e774*/
2794 	         goto TERMINATE;
2795 	
2796 	      /* add facet to the cut and update the rhs and activity of the cut */
2797 	      SCIP_CALL( addFacetToCut(scip, sol, cut, bestfacet, ecaggr->vars, ecaggr->nvars, &cutconstant, &cutactivity,
2798 	            &success) );
2799 	
2800 	      if( !success )
2801 	         goto TERMINATE;
2802 	   }
2803 	
2804 	   /* compute an underestimate for each bilinear term which is not in any edge-concave aggregation */
2805 	   for( i = 0; i < nlrowaggr->nremterms; ++i )
2806 	   {
2807 	      SCIP_VAR* x;
2808 	      SCIP_VAR* y;
2809 	      SCIP_Bool success;
2810 	
2811 	      x = nlrowaggr->remtermvars1[i];
2812 	      y = nlrowaggr->remtermvars2[i];
2813 	      assert(x != NULL);
2814 	      assert(y != NULL);
2815 	
2816 	      SCIP_CALL( addBilinearTermToCut(scip, sol, cut, x, y, nlrowaggr->remtermcoefs[i], &cutconstant, &cutactivity,
2817 	            &success) );
2818 	
2819 	      if( !success )
2820 	         goto TERMINATE;
2821 	   }
2822 	
2823 	   /* add all linear terms to the cut */
2824 	   for( i = 0; i < nlrowaggr->nlinvars; ++i )
2825 	   {
2826 	      SCIP_VAR* x;
2827 	      SCIP_Real coef;
2828 	      SCIP_Bool success;
2829 	
2830 	      x = nlrowaggr->linvars[i];
2831 	      assert(x != NULL);
2832 	
2833 	      coef = nlrowaggr->lincoefs[i];
2834 	
2835 	      SCIP_CALL( addLinearTermToCut(scip, sol, cut, x, coef, &cutconstant, &cutactivity, &success) );
2836 	
2837 	      if( !success )
2838 	         goto TERMINATE;
2839 	   }
2840 	
2841 	   SCIPdebugMsg(scip, "cut activity = %e rhs(nlrow) = %e\n", cutactivity, nlrowaggr->rhs);
2842 	
2843 	   /* set rhs of the cut (substract the constant part of the cut) */
2844 	   SCIP_CALL( SCIPchgRowRhs(scip, cut, nlrowaggr->rhs - cutconstant) );
2845 	   SCIP_CALL( SCIPflushRowExtensions(scip, cut) );
2846 	
2847 	   /* check activity of the row; this assert can fail because of numerics */
2848 	   /* assert(SCIPisFeasEQ(scip, cutactivity - cutconstant, SCIPgetRowSolActivity(scip, cut, sol)) ); */
2849 	
2850 	#ifdef SCIP_DEBUG
2851 	   SCIP_CALL( SCIPprintRow(scip, cut, NULL) );
2852 	#endif
2853 	
2854 	   SCIPdebugMsg(scip, "EC cut <%s>: act=%f eff=%f rank=%d range=%e\n",
2855 	      SCIProwGetName(cut), SCIPgetRowSolActivity(scip, cut, sol), SCIPgetCutEfficacy(scip, sol, cut),
2856 	      SCIProwGetRank(cut), SCIPgetRowMaxCoef(scip, cut) / SCIPgetRowMinCoef(scip, cut) );
2857 	
2858 	   /* try to add the cut has a finite rhs, is efficacious, and does not exceed the maximum cut range */
2859 	   if( !SCIPisInfinity(scip, nlrowaggr->rhs - cutconstant) && SCIPisCutEfficacious(scip, sol, cut)
2860 	      && SCIPgetRowMaxCoef(scip, cut) / SCIPgetRowMinCoef(scip, cut) < sepadata->cutmaxrange )
2861 	   {
2862 	      /* add the cut if it is separating the given solution by at least minviolation */
2863 	      if( SCIPisGE(scip, cutactivity - nlrowaggr->rhs, sepadata->minviolation) )
2864 	      {
2865 	         SCIP_CALL( SCIPaddRow(scip, cut, FALSE, cutoff) );
2866 	         *separated = TRUE;
2867 	         SCIPdebugMsg(scip, "added separating cut\n");
2868 	      }
2869 	
2870 	      if( !(*cutoff) && !islocalcut )
2871 	      {
2872 	         SCIP_CALL( SCIPaddPoolCut(scip, cut) );
2873 	         SCIPdebugMsg(scip, "added cut to cut pool\n");
2874 	      }
2875 	   }
2876 	
2877 	TERMINATE:
2878 	   /* free allocated memory */
2879 	   SCIPfreeBufferArray(scip, &bestfacet);
2880 	
2881 	   /* release the row */
2882 	   SCIP_CALL( SCIPreleaseRow(scip, &cut) );
2883 	
2884 	   return SCIP_OKAY;
2885 	}
2886 	
2887 	/** returns whether it is possible to compute a cut for a given nonlinear row aggregation */
2888 	static
2889 	SCIP_Bool isPossibleToComputeCut(
2890 	   SCIP*                 scip,               /**< SCIP data structure */
2891 	   SCIP_SOL*             sol,                /**< current solution (might be NULL) */
2892 	   SCIP_NLROWAGGR*       nlrowaggr           /**< nonlinear row aggregation */
2893 	   )
2894 	{
2895 	   int i;
2896 	
2897 	   assert(scip != NULL);
2898 	   assert(nlrowaggr != NULL);
2899 	
2900 	   if( !SCIPnlrowIsInNLP(nlrowaggr->nlrow) )
2901 	   {
2902 	      SCIPdebugMsg(scip, "nlrow is not in NLP anymore\n");
2903 	      return FALSE;
2904 	   }
2905 	
2906 	   for( i = 0; i < nlrowaggr->nquadvars; ++i )
2907 	   {
2908 	      SCIP_VAR* var = nlrowaggr->quadvars[i];
2909 	      assert(var != NULL);
2910 	
2911 	      /* check whether the variable has infinite bounds */
2912 	      if( SCIPisInfinity(scip, REALABS(SCIPvarGetLbLocal(var))) || SCIPisInfinity(scip, REALABS(SCIPvarGetUbLocal(var)))
2913 	            || SCIPisInfinity(scip, REALABS(SCIPgetSolVal(scip, sol, var))) )
2914 	      {
2915 	         SCIPdebugMsg(scip, "nlrow aggregation contains unbounded variables\n");
2916 	         return FALSE;
2917 	      }
2918 	
2919 	      /* check whether the variable has been fixed and is in one edge-concave aggregation */
2920 	      if( nlrowaggr->quadvar2aggr[i] >= 0 && SCIPisFeasEQ(scip, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var)) )
2921 	      {
2922 	         SCIPdebugMsg(scip, "nlrow aggregation contains fixed variables in an e.c. aggregation\n");
2923 	         return FALSE;
2924 	      }
2925 	   }
2926 	
2927 	   return TRUE;
2928 	}
2929 	
2930 	/** searches and tries to add edge-concave cuts */
2931 	static
2932 	SCIP_RETCODE separateCuts(
2933 	   SCIP*                 scip,               /**< SCIP data structure */
2934 	   SCIP_SEPA*            sepa,               /**< separator */
2935 	   SCIP_SEPADATA*        sepadata,           /**< separator data */
2936 	   int                   depth,              /**< current depth */
2937 	   SCIP_SOL*             sol,                /**< current solution */
2938 	   SCIP_RESULT*          result              /**< pointer to store the result of the separation call */
2939 	   )
2940 	{
2941 	   int nmaxcuts;
2942 	   int ncuts;
2943 	   int i;
2944 	
2945 	   assert(*result == SCIP_DIDNOTRUN);
2946 	
2947 	   SCIPdebugMsg(scip, "separate cuts...\n");
2948 	
2949 	   /* skip if there are no nonlinear row aggregations */
2950 	   if( sepadata->nnlrowaggrs == 0 )
2951 	   {
2952 	      SCIPdebugMsg(scip, "no aggregations exists -> skip call\n");
2953 	      return SCIP_OKAY;
2954 	   }
2955 	
2956 	   /* get the maximal number of cuts allowed in a separation round */
2957 	   nmaxcuts = depth == 0 ? sepadata->maxsepacutsroot : sepadata->maxsepacuts;
2958 	   ncuts = 0;
2959 	
2960 	   /* try to compute cuts for each nonlinear row independently */
2961 	   for( i = 0; i < sepadata->nnlrowaggrs && ncuts < nmaxcuts && !SCIPisStopped(scip); ++i )
2962 	   {
2963 	      SCIP_NLROWAGGR* nlrowaggr;
2964 	      SCIP_Bool separated;
2965 	      SCIP_Bool cutoff;
2966 	
2967 	      nlrowaggr = sepadata->nlrowaggrs[i];
2968 	      assert(nlrowaggr != NULL);
2969 	
2970 	      /* skip nonlinear aggregations for which it is obviously not possible to compute a cut */
2971 	      if( !isPossibleToComputeCut(scip, sol, nlrowaggr) )
2972 	         return SCIP_OKAY;
2973 	
2974 	      *result = (*result == SCIP_DIDNOTRUN) ? SCIP_DIDNOTFIND : *result;
2975 	
2976 	      SCIPdebugMsg(scip, "try to compute a cut for nonlinear row aggregation %d\n", i);
2977 	
2978 	      /* compute and add cut */
2979 	      SCIP_CALL( computeCut(scip, sepa, sepadata, nlrowaggr, sol, &separated, &cutoff) );
2980 	      SCIPdebugMsg(scip, "found a cut: %s cutoff: %s\n", separated ? "yes" : "no", cutoff ? "yes" : "no");
2981 	
2982 	      /* stop if the current node gets cut off */
2983 	      if( cutoff )
2984 	      {
2985 	         assert(separated);
2986 	         *result = SCIP_CUTOFF;
2987 	         return SCIP_OKAY;
2988 	      }
2989 	
2990 	      /* do not compute more cuts if we already separated the given solution */
2991 	      if( separated )
2992 	      {
2993 	         assert(!cutoff);
2994 	         *result = SCIP_SEPARATED;
2995 	         ++ncuts;
2996 	      }
2997 	   }
2998 	
2999 	   return SCIP_OKAY;
3000 	}
3001 	
3002 	/*
3003 	 * Callback methods of separator
3004 	 */
3005 	
3006 	/** copy method for separator plugins (called when SCIP copies plugins) */
3007 	static
3008 	SCIP_DECL_SEPACOPY(sepaCopyEccuts)
3009 	{  /*lint --e{715}*/
3010 	   assert(scip != NULL);
3011 	   assert(sepa != NULL);
3012 	   assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
3013 	
3014 	   /* call inclusion method of constraint handler */
3015 	   SCIP_CALL( SCIPincludeSepaEccuts(scip) );
3016 	
3017 	   return SCIP_OKAY;
3018 	}
3019 	
3020 	/** destructor of separator to free user data (called when SCIP is exiting) */
3021 	static
3022 	SCIP_DECL_SEPAFREE(sepaFreeEccuts)
3023 	{  /*lint --e{715}*/
3024 	   SCIP_SEPADATA* sepadata;
3025 	
3026 	   sepadata = SCIPsepaGetData(sepa);
3027 	   assert(sepadata != NULL);
3028 	
3029 	   SCIP_CALL( sepadataFree(scip, &sepadata) );
3030 	   SCIPsepaSetData(sepa, NULL);
3031 	
3032 	   return SCIP_OKAY;
3033 	}
3034 	
3035 	/** solving process deinitialization method of separator (called before branch and bound process data is freed) */
3036 	static
3037 	SCIP_DECL_SEPAEXITSOL(sepaExitsolEccuts)
3038 	{  /*lint --e{715}*/
3039 	   SCIP_SEPADATA* sepadata;
3040 	
3041 	   sepadata = SCIPsepaGetData(sepa);
3042 	   assert(sepadata != NULL);
3043 	
3044 	   /* print statistics */
3045 	#ifdef SCIP_STATISTIC
3046 	   SCIPstatisticMessage("rhs-AGGR %d\n", sepadata->nrhsnlrowaggrs);
3047 	   SCIPstatisticMessage("lhs-AGGR %d\n", sepadata->nlhsnlrowaggrs);
3048 	   SCIPstatisticMessage("aggr. search time = %f\n", sepadata->aggrsearchtime);
3049 	#endif
3050 	
3051 	   /* free nonlinear row aggregations */
3052 	   SCIP_CALL( sepadataFreeNlrows(scip, sepadata) );
3053 	
3054 	   /* mark that we should search again for nonlinear row aggregations */
3055 	   sepadata->searchedforaggr = FALSE;
3056 	
3057 	   SCIPdebugMsg(scip, "exitsol\n");
3058 	
3059 	   return SCIP_OKAY;
3060 	}
3061 	
3062 	/** LP solution separation method of separator */
3063 	static
3064 	SCIP_DECL_SEPAEXECLP(sepaExeclpEccuts)
3065 	{  /*lint --e{715}*/
3066 	   SCIP_SEPADATA* sepadata;
3067 	   int ncalls;
3068 	
3069 	   sepadata = SCIPsepaGetData(sepa);
3070 	   assert(sepadata != NULL);
3071 	
3072 	   *result = SCIP_DIDNOTRUN;
3073 	
3074 	   if( !allowlocal )
3075 	      return SCIP_OKAY;
3076 	
3077 	   /* check min- and maximal aggregation size */
3078 	   if( sepadata->maxaggrsize < sepadata->minaggrsize )
3079 	      return SCIP_PARAMETERWRONGVAL;
3080 	
3081 	   /* only call separator, if we are not close to terminating */
3082 	   if( SCIPisStopped(scip) )
3083 	      return SCIP_OKAY;
3084 	
3085 	   /* skip if the LP is not constructed yet */
3086 	   if( !SCIPisNLPConstructed(scip) )
3087 	   {
3088 	      SCIPdebugMsg(scip, "Skip since NLP is not constructed yet.\n");
3089 	      return SCIP_OKAY;
3090 	   }
3091 	
3092 	   /* only call separator up to a maximum depth */
3093 	   if ( sepadata->maxdepth >= 0 && depth > sepadata->maxdepth )
3094 	      return SCIP_OKAY;
3095 	
3096 	   /* only call separator a given number of times at each node */
3097 	   ncalls = SCIPsepaGetNCallsAtNode(sepa);
3098 	   if ( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
3099 	      || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
3100 	      return SCIP_OKAY;
3101 	
3102 	   /* search for nonlinear row aggregations */
3103 	   if( !sepadata->searchedforaggr )
3104 	   {
3105 	      int i;
3106 	
3107 	      SCIPstatistic( sepadata->aggrsearchtime -= SCIPgetTotalTime(scip) );
3108 	
3109 	      SCIPdebugMsg(scip, "search for nonlinear row aggregations\n");
3110 	      for( i = 0; i < SCIPgetNNLPNlRows(scip) && !SCIPisStopped(scip); ++i )
3111 	      {
3112 	         SCIP_NLROW* nlrow = SCIPgetNLPNlRows(scip)[i];
3113 	         SCIP_CALL( findAndStoreEcAggregations(scip, sepadata, nlrow, NULL) );
3114 	      }
3115 	      sepadata->searchedforaggr = TRUE;
3116 	
3117 	      SCIPstatistic( sepadata->aggrsearchtime += SCIPgetTotalTime(scip) );
3118 	   }
3119 	
3120 	   /* search for edge-concave cuts */
3121 	   SCIP_CALL( separateCuts(scip, sepa, sepadata, depth, NULL, result) );
3122 	
3123 	   return SCIP_OKAY;
3124 	}
3125 	
3126 	/*
3127 	 * separator specific interface methods
3128 	 */
3129 	
3130 	/** creates the edge-concave separator and includes it in SCIP
3131 	 *
3132 	 * @ingroup SeparatorIncludes
3133 	 */
3134 	SCIP_RETCODE SCIPincludeSepaEccuts(
3135 	   SCIP*                 scip                /**< SCIP data structure */
3136 	   )
3137 	{
3138 	   SCIP_SEPADATA* sepadata;
3139 	   SCIP_SEPA* sepa;
3140 	
3141 	   /* create eccuts separator data */
3142 	   SCIP_CALL( sepadataCreate(scip, &sepadata) );
3143 	
3144 	   /* include separator */
3145 	   SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
3146 	         SEPA_USESSUBSCIP, SEPA_DELAY, sepaExeclpEccuts, NULL, sepadata) );
3147 	
3148 	   assert(sepa != NULL);
3149 	
3150 	   /* set non fundamental callbacks via setter functions */
3151 	   SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyEccuts) );
3152 	   SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeEccuts) );
3153 	   SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolEccuts) );
3154 	
3155 	   /* add eccuts separator parameters */
3156 	   SCIP_CALL( SCIPaddBoolParam(scip,
3157 	         "separating/" SEPA_NAME "/dynamiccuts",
3158 	         "should generated cuts be removed from the LP if they are no longer tight?",
3159 	         &sepadata->dynamiccuts, FALSE, DEFAULT_DYNAMICCUTS, NULL, NULL) );
3160 	
3161 	   SCIP_CALL( SCIPaddIntParam(scip,
3162 	         "separating/" SEPA_NAME "/maxrounds",
3163 	         "maximal number of eccuts separation rounds per node (-1: unlimited)",
3164 	         &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
3165 	
3166 	   SCIP_CALL( SCIPaddIntParam(scip,
3167 	         "separating/" SEPA_NAME "/maxroundsroot",
3168 	         "maximal number of eccuts separation rounds in the root node (-1: unlimited)",
3169 	         &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
3170 	
3171 	   SCIP_CALL( SCIPaddIntParam(scip,
3172 	         "separating/" SEPA_NAME "/maxdepth",
3173 	         "maximal depth at which the separator is applied (-1: unlimited)",
3174 	         &sepadata->maxdepth, FALSE, DEFAULT_MAXDEPTH, -1, INT_MAX, NULL, NULL) );
3175 	
3176 	   SCIP_CALL( SCIPaddIntParam(scip,
3177 	         "separating/" SEPA_NAME "/maxsepacuts",
3178 	         "maximal number of edge-concave cuts separated per separation round",
3179 	         &sepadata->maxsepacuts, FALSE, DEFAULT_MAXSEPACUTS, 0, INT_MAX, NULL, NULL) );
3180 	
3181 	   SCIP_CALL( SCIPaddIntParam(scip,
3182 	         "separating/" SEPA_NAME "/maxsepacutsroot",
3183 	         "maximal number of edge-concave cuts separated per separation round in the root node",
3184 	         &sepadata->maxsepacutsroot, FALSE, DEFAULT_MAXSEPACUTSROOT, 0, INT_MAX, NULL, NULL) );
3185 	
3186 	   SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/cutmaxrange",
3187 	         "maximal coef. range of a cut (max coef. divided by min coef.) in order to be added to LP relaxation",
3188 	         &sepadata->cutmaxrange, FALSE, DEFAULT_CUTMAXRANGE, 0.0, SCIPinfinity(scip), NULL, NULL) );
3189 	
3190 	   SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/minviolation",
3191 	         "minimal violation of an edge-concave cut to be separated",
3192 	         &sepadata->minviolation, FALSE, DEFAULT_MINVIOLATION, 0.0, 0.5, NULL, NULL) );
3193 	
3194 	   SCIP_CALL( SCIPaddIntParam(scip,
3195 	         "separating/" SEPA_NAME "/minaggrsize",
3196 	         "search for edge-concave aggregations of at least this size",
3197 	         &sepadata->minaggrsize, TRUE, DEFAULT_MINAGGRSIZE, 3, 5, NULL, NULL) );
3198 	
3199 	   SCIP_CALL( SCIPaddIntParam(scip,
3200 	         "separating/" SEPA_NAME "/maxaggrsize",
3201 	         "search for edge-concave aggregations of at most this size",
3202 	         &sepadata->maxaggrsize, TRUE, DEFAULT_MAXAGGRSIZE, 3, 5, NULL, NULL) );
3203 	
3204 	   SCIP_CALL( SCIPaddIntParam(scip,
3205 	         "separating/" SEPA_NAME "/maxbilinterms",
3206 	         "maximum number of bilinear terms allowed to be in a quadratic constraint",
3207 	         &sepadata->maxbilinterms, TRUE, DEFAULT_MAXBILINTERMS, 0, INT_MAX, NULL, NULL) );
3208 	
3209 	   SCIP_CALL( SCIPaddIntParam(scip,
3210 	         "separating/" SEPA_NAME "/maxstallrounds",
3211 	         "maximum number of unsuccessful rounds in the edge-concave aggregation search",
3212 	         &sepadata->maxstallrounds, TRUE, DEFAULT_MAXSTALLROUNDS, 0, INT_MAX, NULL, NULL) );
3213 	
3214 	   return SCIP_OKAY;
3215 	}
3216