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   nlhdlr_soc.c
26   	 * @ingroup DEFPLUGINS_NLHDLR
27   	 * @brief  nonlinear handler for second order cone constraints
28   	
29   	 * @author Benjamin Mueller
30   	 * @author Felipe Serrano
31   	 * @author Fabian Wegscheider
32   	 *
33   	 * This is a nonlinear handler for second order cone constraints of the form
34   	 *
35   	 * \f[\sqrt{\sum_{i=1}^{n} (v_i^T x + \beta_i)^2} \leq v_{n+1}^T x + \beta_{n+1}.\f]
36   	 *
37   	 * Note that \f$v_i\f$, for \f$i \leq n\f$, could be 0, thus allowing a positive constant term inside the root.
38   	 *
39   	 * @todo test if it makes sense to only disaggregate when nterms > some parameter
40   	 *
41   	 */
42   	
43   	#include <string.h>
44   	
45   	#include "scip/nlhdlr_soc.h"
46   	#include "scip/cons_nonlinear.h"
47   	#include "scip/expr_pow.h"
48   	#include "scip/expr_sum.h"
49   	#include "scip/expr_var.h"
50   	#include "scip/debug.h"
51   	#include "scip/pub_nlhdlr.h"
52   	#include "scip/lapack_calls.h"
53   	
54   	
55   	/* fundamental nonlinear handler properties */
56   	#define NLHDLR_NAME               "soc"
57   	#define NLHDLR_DESC               "nonlinear handler for second-order cone structures"
58   	#define NLHDLR_DETECTPRIORITY       100 /**< priority of the nonlinear handler for detection */
59   	#define NLHDLR_ENFOPRIORITY         100 /**< priority of the nonlinear handler for enforcement */
60   	#define DEFAULT_MINCUTEFFICACY     1e-5 /**< default value for parameter mincutefficacy */
61   	#define DEFAULT_COMPEIGENVALUES    TRUE /**< default value for parameter compeigenvalues */
62   	
63   	/*
64   	 * Data structures
65   	 */
66   	
67   	/** nonlinear handler expression data. The data is structured in the following way:
68   	 *
69   	 *  A 'term' is one of the arguments of the quadratic terms, i.e. \f$v_i^T x + beta_i\f$.
70   	 *  The last term is always the one on the right-hand side. This means that nterms is
71   	 *  equal to n+1 in the above description.
72   	 *
73   	 *  - vars contains a list of all expressions which are treated as variables (no duplicates)
74   	 *  - offsets contains the constants beta_i of each term
75   	 *  - transcoefs contains the non-zero values of the transformation vectors v_i of each term
76   	 *  - transcoefsidx contains for each entry of transcoefs the position of the respective variable in vars
77   	 *  - termbegins contains the index at which the transcoefs of each term start, with a sentinel value
78   	 *  - nterms is the total number of terms appearing on both sides
79   	 *  - nvars is the total number of unique variables appearing (length of vars)
80   	 *
81   	 *  Note that the numbers of nonzeroes in v_i is termbegins[i+1] - termbegins[i] and that
82   	 *  the total number of entries in transcoefs and transcoefsidx is termbegins[nterms]
83   	 *
84   	 *  The disaggregation is implicitly stored in the variables disvars and disrow. An SOC as
85   	 *  described above is replaced by n smaller SOCs
86   	 *
87   	 *              (v_i^T x + beta_i)^2 <= disvar_i     * (v_{n+1}^T x + beta_{n+1})
88   	 *
89   	 *  and the row       sum_i disvar_i <= v_{n+1}^T x + beta_{n+1}.
90   	 *
91   	 *  The disaggregation only happens if we have more than 3 terms.
92   	 *
93   	 *  Example: The constraint SQRT(5 + (3x - 4y + 2)^2 + y^2 + 7z^2) <= 5x - y - 1
94   	 *           results in the following nlhdlrexprdata:
95   	 *
96   	 *           vars = {x, y, z}
97   	 *           offsets = {2, 0, 0, SQRT(5), -1}
98   	 *           transcoefs = {3, -4, 1, SQRT(7), 5, -1}
99   	 *           transcoefsidx = {0, 1, 1, 2, 0, 1}
100  	 *           termbegins = {0, 2, 3, 4, 4, 6}
101  	 *           nvars = 3
102  	 *           nterms = 5
103  	 *
104  	 * @note: due to the current implementation, the constant term is the second to last term, except when the SOC was a rotated
105  	 * SOC, e.g., 1 + x^2 - y*z, i.e., when detected by detectSocQuadraticSimple. In that case, the constant is third to
106  	 * last term.
107  	 */
108  	struct SCIP_NlhdlrExprData
109  	{
110  	   SCIP_EXPR**           vars;               /**< expressions which (aux)variables appear on both sides (x) */
111  	   SCIP_Real*            offsets;            /**< offsets of both sides (beta_i) */
112  	   SCIP_Real*            transcoefs;         /**< non-zeros of linear transformation vectors (v_i) */
113  	   int*                  transcoefsidx;      /**< mapping of transformation coefficients to variable indices in vars */
114  	   int*                  termbegins;         /**< starting indices of transcoefs for each term */
115  	   int                   nvars;              /**< total number of variables appearing */
116  	   int                   nterms;             /**< number of summands in the SQRT +1 for RHS (n+1) */
117  	
118  	   /* variables for cone disaggregation */
119  	   SCIP_VAR**            disvars;            /**< disaggregation variables for each term in lhs */
120  	   SCIP_ROW*             disrow;             /**< disaggregation row */
121  	
122  	   /* separation data */
123  	   SCIP_Real*            varvals;            /**< current values for vars */
124  	   SCIP_Real*            disvarvals;         /**< current values for disvars */
125  	};
126  	
127  	struct SCIP_NlhdlrData
128  	{
129  	   SCIP_Real             mincutefficacy;     /**< minimum efficacy a cut need to be added */
130  	   SCIP_Bool             compeigenvalues;    /**< whether Eigenvalue computations should be done to detect complex cases */
131  	};
132  	
133  	/*
134  	 * Local methods
135  	 */
136  	
137  	#ifdef SCIP_DEBUG
138  	/** prints the nlhdlr expression data */
139  	static
140  	void printNlhdlrExprData(
141  	   SCIP*                 scip,               /**< SCIP data structure */
142  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata      /**< pointer to store nonlinear handler expression data */
143  	   )
144  	{
145  	   int nterms;
146  	   int i;
147  	   int j;
148  	
149  	   nterms = nlhdlrexprdata->nterms;
150  	
151  	   SCIPinfoMessage(scip, NULL, "SQRT( ");
152  	
153  	   for( i = 0; i < nterms - 1; ++i )
154  	   {
155  	      int startidx;
156  	
157  	      startidx = nlhdlrexprdata->termbegins[i];
158  	
159  	      /* v_i is 0 */
160  	      if( startidx == nlhdlrexprdata->termbegins[i + 1] )
161  	      {
162  	         assert(nlhdlrexprdata->offsets[i] != 0.0);
163  	
164  	         SCIPinfoMessage(scip, NULL, "%g", SQR(nlhdlrexprdata->offsets[i]));
165  	         continue;
166  	      }
167  	
168  	      /* v_i is not 0 */
169  	      SCIPinfoMessage(scip, NULL, "(");
170  	
171  	      for( j = startidx; j < nlhdlrexprdata->termbegins[i + 1]; ++j )
172  	      {
173  	         if( nlhdlrexprdata->transcoefs[j] != 1.0 )
174  	            SCIPinfoMessage(scip, NULL, " %+g*", nlhdlrexprdata->transcoefs[j]);
175  	         else
176  	            SCIPinfoMessage(scip, NULL, " +");
177  	         if( SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]) != NULL )
178  	         {
179  	            SCIPinfoMessage(scip, NULL, "%s", SCIPvarGetName(SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]])));
180  	            SCIPinfoMessage(scip, NULL, "(%p)", (void*)nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]);
181  	         }
182  	         else
183  	            SCIPinfoMessage(scip, NULL, "%p", (void*)nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]);
184  	      }
185  	      if( nlhdlrexprdata->offsets[i] != 0.0 )
186  	         SCIPinfoMessage(scip, NULL, " %+g", nlhdlrexprdata->offsets[i]);
187  	
188  	      SCIPinfoMessage(scip, NULL, ")^2");
189  	
190  	      if( i < nterms - 2 )
191  	         SCIPinfoMessage(scip, NULL, " + ");
192  	   }
193  	
194  	   SCIPinfoMessage(scip, NULL, " ) <=");
195  	
196  	   for( j = nlhdlrexprdata->termbegins[nterms-1]; j < nlhdlrexprdata->termbegins[nterms]; ++j )
197  	   {
198  	      if( nlhdlrexprdata->transcoefs[j] != 1.0 )
199  	         SCIPinfoMessage(scip, NULL, " %+g*", nlhdlrexprdata->transcoefs[j]);
200  	      else
201  	         SCIPinfoMessage(scip, NULL, " +");
202  	      if( SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]) != NULL )
203  	         SCIPinfoMessage(scip, NULL, "%s", SCIPvarGetName(SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]])));
204  	      else
205  	         SCIPinfoMessage(scip, NULL, "%p", (void*)nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]);
206  	   }
207  	   if( nlhdlrexprdata->offsets[nterms-1] != 0.0 )
208  	      SCIPinfoMessage(scip, NULL, " %+g", nlhdlrexprdata->offsets[nterms-1]);
209  	
210  	   SCIPinfoMessage(scip, NULL, "\n");
211  	}
212  	#endif
213  	
214  	/** helper method to create variables for the cone disaggregation */
215  	static
216  	SCIP_RETCODE createDisaggrVars(
217  	   SCIP*                 scip,               /**< SCIP data structure */
218  	   SCIP_EXPR*            expr,               /**< expression */
219  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata      /**< nonlinear handler expression data */
220  	   )
221  	{
222  	   char name[SCIP_MAXSTRLEN];
223  	   int ndisvars;
224  	   int i;
225  	
226  	   assert(nlhdlrexprdata != NULL);
227  	
228  	   ndisvars = nlhdlrexprdata->nterms - 1;
229  	
230  	   /* allocate memory */
231  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &nlhdlrexprdata->disvars, ndisvars) );
232  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &nlhdlrexprdata->disvarvals, ndisvars) );
233  	
234  	   /* create disaggregation variables representing the epigraph of (v_i^T x + beta_i)^2 / (v_{n+1}^T x + beta_{n+1}) */
235  	   for( i = 0; i < ndisvars; ++i )
236  	   {
237  	      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedis_%p_%d", (void*) expr, i);
238  	      SCIP_CALL( SCIPcreateVarBasic(scip, &nlhdlrexprdata->disvars[i], name, 0.0, SCIPinfinity(scip), 0.0,
239  	               SCIP_VARTYPE_CONTINUOUS) );
240  	      SCIPvarMarkRelaxationOnly(nlhdlrexprdata->disvars[i]);
241  	
242  	      SCIP_CALL( SCIPaddVar(scip, nlhdlrexprdata->disvars[i]) );
243  	      SCIP_CALL( SCIPaddVarLocksType(scip, nlhdlrexprdata->disvars[i], SCIP_LOCKTYPE_MODEL, 1, 1) );
244  	   }
245  	
246  	   return SCIP_OKAY;
247  	}
248  	
249  	/** helper method to free variables for the cone disaggregation */
250  	static
251  	SCIP_RETCODE freeDisaggrVars(
252  	   SCIP*                 scip,               /**< SCIP data structure */
253  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata      /**< nonlinear handler expression data */
254  	   )
255  	{
256  	   int ndisvars;
257  	   int i;
258  	
259  	   assert(nlhdlrexprdata != NULL);
260  	
261  	   if( nlhdlrexprdata->disvars == NULL )
262  	      return SCIP_OKAY;
263  	
264  	   ndisvars = nlhdlrexprdata->nterms - 1;
265  	
266  	   /* release variables */
267  	   for( i = 0; i < ndisvars; ++i )
268  	   {
269  	      SCIP_CALL( SCIPaddVarLocksType(scip, nlhdlrexprdata->disvars[i], SCIP_LOCKTYPE_MODEL, -1, -1) );
270  	      SCIP_CALL( SCIPreleaseVar(scip, &nlhdlrexprdata->disvars[i]) );
271  	   }
272  	
273  	   /* free memory */
274  	   SCIPfreeBlockMemoryArray(scip, &nlhdlrexprdata->disvars, ndisvars);
275  	   SCIPfreeBlockMemoryArrayNull(scip, &nlhdlrexprdata->disvarvals, ndisvars);
276  	
277  	   return SCIP_OKAY;
278  	}
279  	
280  	/** helper method to create the disaggregation row \f$\text{disvars}_i \leq v_{n+1}^T x + \beta_{n+1}\f$ */
281  	static
282  	SCIP_RETCODE createDisaggrRow(
283  	   SCIP*                 scip,               /**< SCIP data structure */
284  	   SCIP_CONSHDLR*        conshdlr,           /**< nonlinear constraint handler */
285  	   SCIP_EXPR*            expr,               /**< expression */
286  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata      /**< nonlinear handler expression data */
287  	   )
288  	{
289  	   SCIP_Real beta;
290  	   char name[SCIP_MAXSTRLEN];
291  	   int ndisvars;
292  	   int nterms;
293  	   int i;
294  	
295  	   assert(scip != NULL);
296  	   assert(expr != NULL);
297  	   assert(nlhdlrexprdata != NULL);
298  	   assert(nlhdlrexprdata->disrow == NULL);
299  	
300  	   nterms = nlhdlrexprdata->nterms;
301  	   beta = nlhdlrexprdata->offsets[nterms - 1];
302  	
303  	   ndisvars = nterms - 1;
304  	
305  	   /* create row 0 <= beta_{n+1} */
306  	   (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedis_%p_row", (void*) expr);
307  	   SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &nlhdlrexprdata->disrow, conshdlr, name,
308  	         -SCIPinfinity(scip), beta, FALSE, FALSE, TRUE) );
309  	
310  	   /* add disvars to row */
311  	   for( i = 0; i < ndisvars; ++i )
312  	   {
313  	      SCIP_CALL( SCIPaddVarToRow(scip, nlhdlrexprdata->disrow, nlhdlrexprdata->disvars[i], 1.0) );
314  	   }
315  	
316  	   /* add rhs vars to row */
317  	   for( i = nlhdlrexprdata->termbegins[nterms - 1]; i < nlhdlrexprdata->termbegins[nterms]; ++i )
318  	   {
319  	      SCIP_VAR* var;
320  	      SCIP_Real coef;
321  	
322  	      var = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[i]]);
323  	      assert(var != NULL);
324  	
325  	      coef = -nlhdlrexprdata->transcoefs[i];
326  	
327  	      SCIP_CALL( SCIPaddVarToRow(scip, nlhdlrexprdata->disrow, var, coef) );
328  	   }
329  	
330  	   return SCIP_OKAY;
331  	}
332  	
333  	/** helper method to create nonlinear handler expression data */
334  	static
335  	SCIP_RETCODE createNlhdlrExprData(
336  	   SCIP*                 scip,               /**< SCIP data structure */
337  	   SCIP_EXPR**           vars,               /**< expressions which variables appear on both sides (\f$x\f$) */
338  	   SCIP_Real*            offsets,            /**< offsets of bot sides (\f$beta_i\f$) */
339  	   SCIP_Real*            transcoefs,         /**< non-zeroes of linear transformation vectors (\f$v_i\f$) */
340  	   int*                  transcoefsidx,      /**< mapping of transformation coefficients to variable indices in vars */
341  	   int*                  termbegins,         /**< starting indices of transcoefs for each term */
342  	   int                   nvars,              /**< total number of variables appearing */
343  	   int                   nterms,             /**< number of summands in the SQRT, +1 for RHS */
344  	   SCIP_NLHDLREXPRDATA** nlhdlrexprdata      /**< pointer to store nonlinear handler expression data */
345  	   )
346  	{
347  	   int ntranscoefs;
348  	
349  	   assert(vars != NULL);
350  	   assert(offsets != NULL);
351  	   assert(transcoefs != NULL);
352  	   assert(transcoefsidx != NULL);
353  	   assert(termbegins != NULL);
354  	   assert(nlhdlrexprdata != NULL);
355  	
356  	   ntranscoefs = termbegins[nterms];
357  	
358  	   SCIP_CALL( SCIPallocBlockMemory(scip, nlhdlrexprdata) );
359  	   SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->vars, vars, nvars) );
360  	   SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->offsets, offsets, nterms) );
361  	   SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->transcoefs, transcoefs, ntranscoefs) );
362  	   SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->transcoefsidx, transcoefsidx, ntranscoefs) );
363  	   SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->termbegins, termbegins, nterms + 1) );
364  	   (*nlhdlrexprdata)->nvars = nvars;
365  	   (*nlhdlrexprdata)->nterms = nterms;
366  	
367  	   (*nlhdlrexprdata)->disrow = NULL;
368  	   (*nlhdlrexprdata)->disvars = NULL;
369  	
370  	   (*nlhdlrexprdata)->varvals = NULL;
371  	   (*nlhdlrexprdata)->disvarvals = NULL;
372  	
373  	#ifdef SCIP_DEBUG
374  	   SCIPdebugMsg(scip, "created nlhdlr data for the following soc expression:\n");
375  	   printNlhdlrExprData(scip, *nlhdlrexprdata);
376  	   /* SCIPdebugMsg(scip, "x is %p\n", (void *)vars[0]); */
377  	#endif
378  	
379  	   return SCIP_OKAY;
380  	}
381  	
382  	/** helper method to free nonlinear handler expression data */
383  	static
384  	SCIP_RETCODE freeNlhdlrExprData(
385  	   SCIP*                 scip,               /**< SCIP data structure */
386  	   SCIP_NLHDLREXPRDATA** nlhdlrexprdata      /**< pointer to free nonlinear handler expression data */
387  	   )
388  	{
389  	   int ntranscoefs;
390  	
391  	   assert(nlhdlrexprdata != NULL);
392  	   assert(*nlhdlrexprdata != NULL);
393  	
394  	   /* free variables and row for cone disaggregation */
395  	   SCIP_CALL( freeDisaggrVars(scip, *nlhdlrexprdata) );
396  	
397  	   ntranscoefs = (*nlhdlrexprdata)->termbegins[(*nlhdlrexprdata)->nterms];
398  	
399  	   SCIPfreeBlockMemoryArrayNull(scip, &(*nlhdlrexprdata)->varvals, (*nlhdlrexprdata)->nvars);
400  	   SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->termbegins, (*nlhdlrexprdata)->nterms + 1);
401  	   SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->transcoefsidx, ntranscoefs);
402  	   SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->transcoefs, ntranscoefs);
403  	   SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->offsets, (*nlhdlrexprdata)->nterms);
404  	   SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->vars, (*nlhdlrexprdata)->nvars);
405  	   SCIPfreeBlockMemory(scip, nlhdlrexprdata);
406  	
407  	   return SCIP_OKAY;
408  	}
409  	
410  	/** set varvalrs in nlhdlrexprdata to values from given SCIP solution */
411  	static
412  	void updateVarVals(
413  	   SCIP*                 scip,               /**< SCIP data structure */
414  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nonlinear handler expression data */
415  	   SCIP_SOL*             sol,                /**< SCIP solution */
416  	   SCIP_Bool             roundtinyfrac       /**< whether values close to integers should be rounded */
417  	   )
418  	{
419  	   int i;
420  	
421  	   assert(nlhdlrexprdata != NULL);
422  	   assert(nlhdlrexprdata->varvals != NULL);
423  	
424  	   /* update varvals */
425  	   for( i = 0; i < nlhdlrexprdata->nvars; ++i )
426  	   {
427  	      nlhdlrexprdata->varvals[i] = SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[i]));
428  	      if( roundtinyfrac && SCIPisIntegral(scip, nlhdlrexprdata->varvals[i]) )
429  	         nlhdlrexprdata->varvals[i] = SCIPround(scip, nlhdlrexprdata->varvals[i]);
430  	   }
431  	
432  	   /* update disvarvals (in unittests, this may be NULL even though nterms > 1 */
433  	   if( nlhdlrexprdata->disvarvals != NULL )
434  	      for( i = 0; i < nlhdlrexprdata->nterms - 1; ++i )
435  	      {
436  	         nlhdlrexprdata->disvarvals[i] = SCIPgetSolVal(scip, sol, nlhdlrexprdata->disvars[i]);
437  	         if( roundtinyfrac && SCIPisIntegral(scip, nlhdlrexprdata->disvarvals[i]) )
438  	            nlhdlrexprdata->disvarvals[i] = SCIPround(scip, nlhdlrexprdata->disvarvals[i]);
439  	      }
440  	}
441  	
442  	/** evaluate a single term of the form \f$v_i^T x + \beta_i\f$ */
443  	static
444  	SCIP_Real evalSingleTerm(
445  	   SCIP*                 scip,               /**< SCIP data structure */
446  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nonlinear handler expression data */
447  	   int                   k                   /**< term to be evaluated */
448  	   )
449  	{
450  	   SCIP_Real result;
451  	   int i;
452  	
453  	   assert(scip != NULL);
454  	   assert(nlhdlrexprdata != NULL);
455  	   assert(0 <= k && k < nlhdlrexprdata->nterms);
456  	
457  	   result = nlhdlrexprdata->offsets[k];
458  	
459  	   for( i = nlhdlrexprdata->termbegins[k]; i < nlhdlrexprdata->termbegins[k + 1]; ++i )
460  	      result += nlhdlrexprdata->transcoefs[i] * nlhdlrexprdata->varvals[nlhdlrexprdata->transcoefsidx[i]];
461  	
462  	   return result;
463  	}
464  	
465  	/** computes gradient cut for a 2D or 3D SOC
466  	 *
467  	 *  A 3D SOC looks like
468  	 *  \f[
469  	 *    \sqrt{ (v_1^T x + \beta_1)^2 + (v_2^T x + \beta_2)^2 } \leq v_3^T x + \beta_3
470  	 *  \f]
471  	 *
472  	 *  Let \f$f(x)\f$ be the left-hand-side. The partial derivatives of \f$f\f$ are given by
473  	 *  \f[
474  	 *    \frac{\delta f}{\delta x_j} = \frac{(v_1)_j(v_1^T x + \beta_1) + (v_2)_j (v_2^T x + \beta_2)}{f(x)}
475  	 *  \f]
476  	 *
477  	 *  and the gradient cut is then \f$f(x^*) + \nabla f(x^*)(x - x^*) \leq v_3^T x + \beta_3\f$.
478  	 *
479  	 *  If \f$\beta_1 = \beta_2 = 0\f$, then the constant on the left-hand-side of the cut becomes zero:
480  	 *  \f[
481  	 *    f(x^*) - (\frac{(v_1)_j v_1^T x^* + (v_2)_j v_2^T x^*}{f(x^*)})_j^T x^*
482  	 *    = f(x^*) - \frac{1}{f(x^*)} \sum_j ((v_1)_j x_j^* v_1^T x^* + (v_2)_j x_j^* v_2^T x^*)
483  	 *    = f(x^*) - \frac{1}{f(x^*)} ((v_1^T x^*)^2 + (v_2^T x^*)^2)
484  	 *    = f(x^*) - \frac{1}{f(x^*)} f(x^*)^2 = 0
485  	 *  \f]
486  	 *
487  	 *  A 2D SOC is
488  	 *  \f[
489  	 *    |v_1^T x + \beta_1| \leq v_2^T x + \beta_2
490  	 *  \f]
491  	 *  but we build the cut using the same procedure as for 3D.
492  	 */
493  	static
494  	SCIP_RETCODE generateCutSolSOC(
495  	   SCIP*                 scip,               /**< SCIP data structure */
496  	   SCIP_ROWPREP**        rowprep,            /**< buffer to store rowprep with cut data */
497  	   SCIP_EXPR*            expr,               /**< expression */
498  	   SCIP_CONS*            cons,               /**< the constraint that expr is part of */
499  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nonlinear handler expression data */
500  	   SCIP_Real             mincutviolation,    /**< minimal required cut violation */
501  	   SCIP_Real             rhsval              /**< value of last term at sol */
502  	   )
503  	{
504  	   SCIP_Real* transcoefs;
505  	   SCIP_Real cutcoef;
506  	   SCIP_Real fvalue;
507  	   SCIP_Real valterms[2] = {0.0, 0.0}; /* for lint */
508  	   SCIP_Real cutrhs;
509  	   SCIP_EXPR** vars;
510  	   SCIP_VAR* cutvar;
511  	   SCIP_Bool offsetzero;
512  	   int* transcoefsidx;
513  	   int* termbegins;
514  	   int nterms;
515  	   int i;
516  	   int j;
517  	
518  	   assert(rowprep != NULL);
519  	   assert(expr != NULL);
520  	   assert(cons != NULL);
521  	   assert(nlhdlrexprdata != NULL);
522  	
523  	   vars = nlhdlrexprdata->vars;
524  	   transcoefs = nlhdlrexprdata->transcoefs;
525  	   transcoefsidx = nlhdlrexprdata->transcoefsidx;
526  	   termbegins = nlhdlrexprdata->termbegins;
527  	   nterms = nlhdlrexprdata->nterms;
528  	
529  	   *rowprep = NULL;
530  	
531  	   /* evaluate lhs terms and compute f(x*), check whether both beta_1 and beta_2 are zero */
532  	   fvalue = 0.0;
533  	   offsetzero = TRUE;
534  	   for( i = 0; i < nterms - 1; ++i )
535  	   {
536  	      valterms[i] = evalSingleTerm(scip, nlhdlrexprdata, i);
537  	      fvalue += SQR( valterms[i] );
538  	      if( nlhdlrexprdata->offsets[i] != 0.0 )
539  	         offsetzero = FALSE;
540  	   }
541  	   fvalue = SQRT( fvalue );
542  	
543  	   /* don't generate cut if we are not violated @todo: remove this once core detects better when a nlhdlr's cons is
544  	    * violated
545  	    */
546  	   if( fvalue - rhsval <= mincutviolation )
547  	   {
548  	      SCIPdebugMsg(scip, "do not generate cut: rhsval %g, fvalue %g violation is %g\n", rhsval, fvalue, fvalue - rhsval);
549  	      return SCIP_OKAY;
550  	   }
551  	
552  	   /* if f(x*) = 0 then we are at top of cone, where we cannot generate cut */
553  	   if( SCIPisZero(scip, fvalue) )
554  	   {
555  	      SCIPdebugMsg(scip, "do not generate cut for lhs=%g, cannot linearize at top of cone\n", fvalue);
556  	      return SCIP_OKAY;
557  	   }
558  	
559  	   /* create cut */
560  	   SCIP_CALL( SCIPcreateRowprep(scip, rowprep, SCIP_SIDETYPE_RIGHT, FALSE) );
561  	   SCIP_CALL( SCIPensureRowprepSize(scip, *rowprep, termbegins[nterms]) );
562  	
563  	   /* cut is f(x*) + \nabla f(x*)^T (x - x*) \leq v_n^T x + \beta_n, i.e.,
564  	    * \nabla f(x*)^T x - v_n^T x \leq \beta_n + \nabla f(x*)^T x* - f(x*)
565  	    * thus cutrhs is \beta_n - f(x*) + \nabla f(x*)^T x*
566  	    * if offsetzero, then we make sure that cutrhs is exactly \beta_n
567  	    */
568  	   cutrhs = nlhdlrexprdata->offsets[nterms - 1];
569  	   if( !offsetzero )
570  	      cutrhs -= fvalue;
571  	
572  	   /* add cut coefficients from lhs terms and compute cut's rhs */
573  	   for( j = 0; j < nterms - 1; ++j )
574  	   {
575  	      for( i = termbegins[j]; i < termbegins[j + 1]; ++i )
576  	      {
577  	         cutvar = SCIPgetExprAuxVarNonlinear(vars[transcoefsidx[i]]);
578  	
579  	         /* cutcoef is (the first part of) the partial derivative w.r.t cutvar */
580  	         cutcoef = transcoefs[i] * valterms[j] / fvalue;
581  	
582  	         SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, cutvar, cutcoef) );
583  	
584  	         if( !offsetzero )
585  	            cutrhs += cutcoef * nlhdlrexprdata->varvals[transcoefsidx[i]];
586  	      }
587  	   }
588  	
589  	   /* add terms for v_n */
590  	   for( i = termbegins[nterms - 1]; i < termbegins[nterms]; ++i )
591  	   {
592  	      cutvar = SCIPgetExprAuxVarNonlinear(vars[transcoefsidx[i]]);
593  	      SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, cutvar, -transcoefs[i]) );
594  	   }
595  	
596  	   /* add side */
597  	   SCIProwprepAddSide(*rowprep, cutrhs);
598  	
599  	   /* set name */
(1) Event invalid_type: Argument "SCIPgetNLPs(scip)" to format specifier "%d" was expected to have type "int" but has type "long long". [details]
600  	   (void) SCIPsnprintf(SCIProwprepGetName(*rowprep), SCIP_MAXSTRLEN, "soc%d_%p_%d", nterms, (void*) expr, SCIPgetNLPs(scip));
601  	
602  	   return SCIP_OKAY;
603  	}
604  	
605  	/** helper method to compute and add a gradient cut for the k-th cone disaggregation
606  	 *
607  	 *  After the SOC constraint \f$\sqrt{\sum_{i = 0}^{n-1} (v_i^T x + \beta_i)^2} \leq v_n^T x + \beta_n\f$
608  	 *  has been disaggregated into the row \f$\sum_{i = 0}^{n-1} y_i \leq v_n^T x + \beta_n\f$ and the smaller SOC constraints
609  	 *  \f[
610  	 *    (v_i^T x + \beta_i)^2 \leq (v_n^T x + \beta_n) y_i \text{ for } i \in \{0, \ldots, n -1\},
611  	 *  \f]
612  	 *  we want to separate one of the small rotated cones.
613  	 *  We first transform it into standard form:
614  	 *  \f[
615  	 *    \sqrt{4(v_i^T x + \beta_i)^2 + (v_n^T x + \beta_n - y_i)^2} - v_n^T x - \beta_n - y_i \leq 0.
616  	 *  \f]
617  	 *  Let \f$f(x,y)\f$ be the left-hand-side. We now compute the gradient by
618  	 *  \f{align*}{
619  	 *    \frac{\delta f}{\delta x_j} &= \frac{(v_i)_j(4v_i^T x + 4\beta_i) + (v_n)_j(v_n^T x + \beta_n - y_i)}{\sqrt{4(v_i^T x + \beta_i)^2 + (v_n^T x + \beta_n - y_i)^2}} - (v_n)_j \\
620  	 *    \frac{\delta f}{\delta y_i} &= \frac{y_i - v_n^T x -\beta_n}{\sqrt{4(v_i^T x + \beta_i)^2 + (v_n^T x + \beta_n - y_i)^2}} - 1
621  	 *  \f}
622  	 *  and the gradient cut is then \f$f(x^*, y^*) + \nabla f(x^*,y^*)((x,y) - (x^*, y^*)) \leq 0\f$.
623  	 *
624  	 *  As in \ref generateCutSolSOC(), the cut constant is zero if \f$\beta_i = \beta_n = 0\f$.
625  	 */
626  	static
627  	SCIP_RETCODE generateCutSolDisagg(
628  	   SCIP*                 scip,               /**< SCIP data structure */
629  	   SCIP_ROWPREP**        rowprep,            /**< buffer to store rowprep with cut data */
630  	   SCIP_EXPR*            expr,               /**< expression */
631  	   SCIP_CONS*            cons,               /**< the constraint that expr is part of */
632  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nonlinear handler expression data */
633  	   int                   disaggidx,          /**< index of disaggregation to separate */
634  	   SCIP_Real             mincutviolation,    /**< minimal required cut violation */
635  	   SCIP_Real             rhsval              /**< value of the rhs term */
636  	   )
637  	{
638  	   SCIP_EXPR** vars;
639  	   SCIP_VAR** disvars;
640  	   SCIP_Real* transcoefs;
641  	   int* transcoefsidx;
642  	   int* termbegins;
643  	   SCIP_VAR* cutvar;
644  	   SCIP_Real cutcoef;
645  	   SCIP_Real fvalue;
646  	   SCIP_Real disvarval;
647  	   SCIP_Real lhsval;
648  	   SCIP_Real constant;
649  	   SCIP_Real denominator;
650  	   SCIP_Bool offsetzero;
651  	   int ncutvars;
652  	   int nterms;
653  	   int i;
654  	
655  	   assert(rowprep != NULL);
656  	   assert(expr != NULL);
657  	   assert(cons != NULL);
658  	   assert(nlhdlrexprdata != NULL);
659  	   assert(disaggidx < nlhdlrexprdata->nterms-1);
660  	
661  	   vars = nlhdlrexprdata->vars;
662  	   disvars = nlhdlrexprdata->disvars;
663  	   transcoefs = nlhdlrexprdata->transcoefs;
664  	   transcoefsidx = nlhdlrexprdata->transcoefsidx;
665  	   termbegins = nlhdlrexprdata->termbegins;
666  	   nterms = nlhdlrexprdata->nterms;
667  	
668  	   /* nterms is equal to n in the description and disaggidx is in {0, ..., n - 1} */
669  	
670  	   *rowprep = NULL;
671  	
672  	   disvarval = nlhdlrexprdata->disvarvals[disaggidx];
673  	
674  	   lhsval = evalSingleTerm(scip, nlhdlrexprdata, disaggidx);
675  	
676  	   denominator = SQRT(4.0 * SQR(lhsval) + SQR(rhsval - disvarval));
677  	
678  	   /* compute value of function to be separated (f(x*,y*)) */
679  	   fvalue = denominator - rhsval - disvarval;
680  	
681  	   /* if the disagg soc is not violated don't compute cut */
682  	   if( fvalue <= mincutviolation )
683  	   {
684  	      SCIPdebugMsg(scip, "skip cut on disaggregation index %d as violation=%g below minviolation %g\n", disaggidx,
685  	            fvalue, mincutviolation);
686  	      return SCIP_OKAY;
687  	   }
688  	
689  	   /* if the denominator is 0 -> the constraint can't be violated, and the gradient is infinite */
690  	   if( SCIPisZero(scip, denominator) )
691  	   {
692  	      SCIPdebugMsg(scip, "skip cut on disaggregation index %d as we are on top of cone (denom=%g)\n", disaggidx, denominator);
693  	      return SCIP_OKAY;
694  	   }
695  	
696  	   /* compute upper bound on the number of variables in cut: vars in rhs + vars in term + disagg var */
697  	   ncutvars = (termbegins[nterms] - termbegins[nterms-1]) + (termbegins[disaggidx + 1] - termbegins[disaggidx]) + 1;
698  	
699  	   /* create cut */
700  	   SCIP_CALL( SCIPcreateRowprep(scip, rowprep, SCIP_SIDETYPE_RIGHT, FALSE) );
701  	   SCIP_CALL( SCIPensureRowprepSize(scip, *rowprep, ncutvars) );
702  	
703  	   /* check whether offsets (beta) are zero, so we can know cut constant will be zero */
704  	   offsetzero = nlhdlrexprdata->offsets[disaggidx] == 0.0 && nlhdlrexprdata->offsets[nterms-1] == 0.0;
705  	
706  	   /* constant will be grad_f(x*,y*)^T  (x*, y*) */
707  	   constant = 0.0;
708  	
709  	   /* a variable could appear on the lhs and rhs, but we add the coefficients separately  */
710  	
711  	   /* add terms for v_disaggidx */
712  	   for( i = termbegins[disaggidx]; i < termbegins[disaggidx + 1]; ++i )
713  	   {
714  	      cutvar = SCIPgetExprAuxVarNonlinear(vars[transcoefsidx[i]]);
715  	      assert(cutvar != NULL);
716  	
717  	      /* cutcoef is (the first part of) the partial derivative w.r.t cutvar */
718  	      cutcoef = 4.0 * lhsval * transcoefs[i] / denominator;
719  	
720  	      SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, cutvar, cutcoef) );
721  	
722  	      if( !offsetzero )
723  	         constant += cutcoef * nlhdlrexprdata->varvals[transcoefsidx[i]];
724  	   }
725  	
726  	   /* add terms for v_n */
727  	   for( i = termbegins[nterms - 1]; i < termbegins[nterms]; ++i )
728  	   {
729  	      cutvar = SCIPgetExprAuxVarNonlinear(vars[transcoefsidx[i]]);
730  	      assert(cutvar != NULL);
731  	
732  	      /* cutcoef is the (second part of) the partial derivative w.r.t cutvar */
733  	      cutcoef = (rhsval - disvarval) * transcoefs[i] / denominator - transcoefs[i];
734  	
735  	      SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, cutvar, cutcoef) );
736  	
737  	      if( !offsetzero )
738  	         constant += cutcoef * nlhdlrexprdata->varvals[transcoefsidx[i]];
739  	   }
740  	
741  	   /* add term for disvar: cutcoef is the the partial derivative w.r.t. the disaggregation variable */
742  	   cutcoef = (disvarval - rhsval) / denominator - 1.0;
743  	   cutvar = disvars[disaggidx];
744  	
745  	   SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, cutvar, cutcoef) );
746  	
747  	   if( !offsetzero )
748  	   {
749  	      constant += cutcoef * nlhdlrexprdata->disvarvals[disaggidx];
750  	
751  	      /* add side */
752  	      SCIProwprepAddSide(*rowprep, constant - fvalue);
753  	   }
754  	
755  	   /* set name */
756  	   (void) SCIPsnprintf(SCIProwprepGetName(*rowprep), SCIP_MAXSTRLEN, "soc_%p_%d_%d", (void*) expr, disaggidx, SCIPgetNLPs(scip));
757  	
758  	   return SCIP_OKAY;
759  	}
760  	
761  	/** given a rowprep, does a number of cleanup and checks and, if successful, generate a cut to be added to the sepastorage */
762  	static
763  	SCIP_RETCODE addCut(
764  	   SCIP*                 scip,               /**< SCIP data structure */
765  	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< nonlinear handler data */
766  	   SCIP_ROWPREP*         rowprep,            /**< rowprep from which to generate row and add as cut */
767  	   SCIP_SOL*             sol,                /**< solution to be separated */
768  	   SCIP_CONS*            cons,               /**< constraint for which cut is generated, or NULL */
769  	   SCIP_Bool             allowweakcuts,      /**< whether weak cuts are allowed */
770  	   SCIP_RESULT*          result              /**< result pointer to update (set to SCIP_CUTOFF or SCIP_SEPARATED if cut is added) */
771  	   )
772  	{
773  	   SCIP_ROW* cut;
774  	   SCIP_Real cutefficacy;
775  	   SCIP_Bool success;
776  	
777  	   assert(scip != NULL);
778  	   assert(nlhdlrdata != NULL);
779  	   assert(rowprep != NULL);
780  	   assert(result != NULL);
781  	
782  	   SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, sol, SCIPgetHugeValue(scip), &success) );
783  	
784  	   if( !success )
785  	   {
786  	      SCIPdebugMsg(scip, "rowprep cleanup failed, skip cut\n");
787  	      return SCIP_OKAY;
788  	   }
789  	
790  	   if( SCIPgetRowprepViolation(scip, rowprep, sol, NULL) <= SCIPgetLPFeastol(scip) )
791  	   {
792  	      SCIPdebugMsg(scip, "rowprep violation %g below LP feastol %g, skip cut\n",
793  	         SCIPgetRowprepViolation(scip, rowprep, sol, NULL), SCIPgetLPFeastol(scip));
794  	      return SCIP_OKAY;
795  	   }
796  	
797  	   SCIP_CALL( SCIPgetRowprepRowCons(scip, &cut, rowprep, cons) );
798  	
799  	   cutefficacy = SCIPgetCutEfficacy(scip, sol, cut);
800  	
801  	   SCIPdebugMsg(scip, "generated row for SOC, efficacy=%g, minefficacy=%g, allowweakcuts=%u\n",
802  	      cutefficacy, nlhdlrdata->mincutefficacy, allowweakcuts);
803  	
804  	   /* check whether cut is applicable */
805  	   if( SCIPisCutApplicable(scip, cut) && (allowweakcuts || cutefficacy >= nlhdlrdata->mincutefficacy) )
806  	   {
807  	      SCIP_Bool infeasible;
808  	
809  	      SCIP_CALL( SCIPaddRow(scip, cut, FALSE, &infeasible) );
810  	
811  	#ifdef SCIP_CONSNONLINEAR_ROWNOTREMOVABLE
812  	      /* mark row as not removable from LP for current node, if in enforcement (==addbranchscores) (this can prevent some cycling) */
813  	      if( addbranchscores )
814  	         SCIPmarkRowNotRemovableLocal(scip, row);
815  	#endif
816  	
817  	      if( infeasible )
818  	         *result = SCIP_CUTOFF;
819  	      else
820  	         *result = SCIP_SEPARATED;
821  	   }
822  	
823  	   /* release row */
824  	   SCIP_CALL( SCIPreleaseRow(scip, &cut) );
825  	
826  	   return SCIP_OKAY;
827  	}
828  	
829  	/** given a rowprep, does a number of cleanup and checks and, if successful, generate a cut to be added to the cutpool */
830  	static
831  	SCIP_RETCODE addCutPool(
832  	   SCIP*                 scip,               /**< SCIP data structure */
833  	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< nonlinear handler data */
834  	   SCIP_ROWPREP*         rowprep,            /**< rowprep from which to generate row and add as cut */
835  	   SCIP_SOL*             sol,                /**< solution to be separated */
836  	   SCIP_CONS*            cons                /**< constraint for which cut is generated, or NULL */
837  	   )
838  	{
839  	   SCIP_ROW* cut;
840  	   SCIP_Bool success;
841  	
842  	   assert(scip != NULL);
843  	   assert(nlhdlrdata != NULL);
844  	   assert(rowprep != NULL);
845  	
846  	   assert(!SCIProwprepIsLocal(rowprep));
847  	
848  	   SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, sol, SCIPgetHugeValue(scip), &success) );
849  	   /* if failed or cut is only locally valid now, then skip */
850  	   if( !success || SCIProwprepIsLocal(rowprep) )
851  	      return SCIP_OKAY;
852  	
853  	   /* if row after cleanup is just a boundchange, then skip */
854  	   if( SCIProwprepGetNVars(rowprep) <= 1 )
855  	      return SCIP_OKAY;
856  	
857  	   /* generate row and add to cutpool */
858  	   SCIP_CALL( SCIPgetRowprepRowCons(scip, &cut, rowprep, cons) );
859  	
860  	   SCIP_CALL( SCIPaddPoolCut(scip, cut) );
861  	
862  	   SCIP_CALL( SCIPreleaseRow(scip, &cut) );
863  	
864  	   return SCIP_OKAY;
865  	}
866  	
867  	/** checks if an expression is quadratic and collects all occurring expressions
868  	 *
869  	 * @pre `expr2idx` and `occurringexprs` need to be initialized with capacity 2 * nchildren
870  	 *
871  	 * @note We assume that a linear term always appears before its corresponding
872  	 * quadratic term in quadexpr; this should be ensured by canonicalize
873  	 */
874  	static
875  	SCIP_RETCODE checkAndCollectQuadratic(
876  	   SCIP*                 scip,               /**< SCIP data structure */
877  	   SCIP_EXPR*            quadexpr,           /**< candidate for a quadratic expression */
878  	   SCIP_HASHMAP*         expr2idx,           /**< hashmap to store expressions */
879  	   SCIP_EXPR**           occurringexprs,     /**< array to store expressions */
880  	   int*                  nexprs,             /**< buffer to store number of expressions */
881  	   SCIP_Bool*            success             /**< buffer to store whether the check was successful */
882  	   )
883  	{
884  	   SCIP_EXPR** children;
885  	   int nchildren;
886  	   int i;
887  	
888  	   assert(scip != NULL);
889  	   assert(quadexpr != NULL);
890  	   assert(expr2idx != NULL);
891  	   assert(occurringexprs != NULL);
892  	   assert(nexprs != NULL);
893  	   assert(success != NULL);
894  	
895  	   *nexprs = 0;
896  	   *success = FALSE;
897  	   children = SCIPexprGetChildren(quadexpr);
898  	   nchildren = SCIPexprGetNChildren(quadexpr);
899  	
900  	   /* iterate in reverse order to ensure that quadratic terms are found before linear terms */
901  	   for( i = nchildren - 1; i >= 0; --i )
902  	   {
903  	      SCIP_EXPR* child;
904  	
905  	      child = children[i];
906  	      if( SCIPisExprPower(scip, child) )
907  	      {
908  	         SCIP_EXPR* childarg;
909  	
910  	         if( SCIPgetExponentExprPow(child) != 2.0 )
911  	            return SCIP_OKAY;
912  	
913  	         childarg = SCIPexprGetChildren(child)[0];
914  	
915  	         if( !SCIPhashmapExists(expr2idx, (void*) childarg) )
916  	         {
917  	            SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void*) childarg, *nexprs) );
918  	
919  	            /* store the expression so we know it later */
920  	            assert(*nexprs < 2 * nchildren);
921  	            occurringexprs[*nexprs] = childarg;
922  	
923  	            ++(*nexprs);
924  	         }
925  	      }
926  	      else if( SCIPisExprVar(scip, child) && SCIPvarIsBinary(SCIPgetVarExprVar(child)) )
927  	      {
928  	         if( !SCIPhashmapExists(expr2idx, (void*) child) )
929  	         {
930  	            SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void*) child, *nexprs) );
931  	
932  	            /* store the expression so we know it later */
933  	            assert(*nexprs < 2 * nchildren);
934  	            occurringexprs[*nexprs] = child;
935  	
936  	            ++(*nexprs);
937  	         }
938  	      }
939  	      else if( SCIPisExprProduct(scip, child) )
940  	      {
941  	         SCIP_EXPR* childarg1;
942  	         SCIP_EXPR* childarg2;
943  	
944  	         if( SCIPexprGetNChildren(child) != 2 )
945  	            return SCIP_OKAY;
946  	
947  	         childarg1 = SCIPexprGetChildren(child)[0];
948  	         childarg2 = SCIPexprGetChildren(child)[1];
949  	
950  	         if( !SCIPhashmapExists(expr2idx, (void*) childarg1) )
951  	         {
952  	            SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void*) childarg1, *nexprs) );
953  	
954  	            /* store the expression so we know it later */
955  	            assert(*nexprs < 2 * nchildren);
956  	            occurringexprs[*nexprs] = childarg1;
957  	
958  	            ++(*nexprs);
959  	         }
960  	
961  	         if( !SCIPhashmapExists(expr2idx, (void*) childarg2) )
962  	         {
963  	            SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void*) childarg2, *nexprs) );
964  	
965  	            /* store the expression so we know it later */
966  	            assert(*nexprs < 2 * nchildren);
967  	            occurringexprs[*nexprs] = childarg2;
968  	
969  	            ++(*nexprs);
970  	         }
971  	      }
972  	      else
973  	      {
974  	         /* if there is a linear term without corresponding quadratic term, it is not a SOC */
975  	         if( !SCIPhashmapExists(expr2idx, (void*) child) )
976  	            return SCIP_OKAY;
977  	      }
978  	   }
979  	
980  	   *success = TRUE;
981  	
982  	   return SCIP_OKAY;
983  	}
984  	
985  	/* builds the constraint defining matrix and vector of a quadratic expression
986  	 *
987  	 * @pre `quadmatrix` and `linvector` need to be initialized with size `nexprs`^2 and `nexprs`, resp.
988  	 */
989  	static
990  	void buildQuadExprMatrix(
991  	   SCIP*                 scip,               /**< SCIP data structure */
992  	   SCIP_EXPR*            quadexpr,           /**< the quadratic expression */
993  	   SCIP_HASHMAP*         expr2idx,           /**< hashmap mapping the occurring expressions to their index */
994  	   int                   nexprs,             /**< number of occurring expressions */
995  	   SCIP_Real*            quadmatrix,         /**< pointer to store (the lower-left triangle of) the quadratic matrix */
996  	   SCIP_Real*            linvector           /**< pointer to store the linear vector */
997  	   )
998  	{
999  	   SCIP_EXPR** children;
1000 	   SCIP_Real* childcoefs;
1001 	   int nchildren;
1002 	   int i;
1003 	
1004 	   assert(scip != NULL);
1005 	   assert(quadexpr != NULL);
1006 	   assert(expr2idx != NULL);
1007 	   assert(quadmatrix != NULL);
1008 	   assert(linvector != NULL);
1009 	
1010 	   children = SCIPexprGetChildren(quadexpr);
1011 	   nchildren = SCIPexprGetNChildren(quadexpr);
1012 	   childcoefs = SCIPgetCoefsExprSum(quadexpr);
1013 	
1014 	   /* iterate over children to build the constraint defining matrix and vector */
1015 	   for( i = 0; i < nchildren; ++i )
1016 	   {
1017 	      int varpos;
1018 	
1019 	      if( SCIPisExprPower(scip, children[i]) )
1020 	      {
1021 	         assert(SCIPgetExponentExprPow(children[i]) == 2.0);
1022 	         assert(SCIPhashmapExists(expr2idx, (void*) SCIPexprGetChildren(children[i])[0]));
1023 	
1024 	         varpos = SCIPhashmapGetImageInt(expr2idx, (void*) SCIPexprGetChildren(children[i])[0]);
1025 	         assert(0 <= varpos && varpos < nexprs);
1026 	
1027 	         quadmatrix[varpos * nexprs + varpos] = childcoefs[i];
1028 	      }
1029 	      else if( SCIPisExprVar(scip, children[i]) && SCIPvarIsBinary(SCIPgetVarExprVar(children[i])) )
1030 	      {
1031 	         assert(SCIPhashmapExists(expr2idx, (void*) children[i]));
1032 	
1033 	         varpos = SCIPhashmapGetImageInt(expr2idx, (void*) children[i]);
1034 	         assert(0 <= varpos && varpos < nexprs);
1035 	
1036 	         quadmatrix[varpos * nexprs + varpos] = childcoefs[i];
1037 	      }
1038 	      else if( SCIPisExprProduct(scip, children[i]) )
1039 	      {
1040 	         int varpos2;
1041 	
1042 	         assert(SCIPexprGetNChildren(children[i]) == 2);
1043 	         assert(SCIPhashmapExists(expr2idx, (void*) SCIPexprGetChildren(children[i])[0]));
1044 	         assert(SCIPhashmapExists(expr2idx, (void*) SCIPexprGetChildren(children[i])[1]));
1045 	
1046 	         varpos = SCIPhashmapGetImageInt(expr2idx, (void*) SCIPexprGetChildren(children[i])[0]);
1047 	         assert(0 <= varpos && varpos < nexprs);
1048 	
1049 	         varpos2 = SCIPhashmapGetImageInt(expr2idx, (void*) SCIPexprGetChildren(children[i])[1]);
1050 	         assert(0 <= varpos2 && varpos2 < nexprs);
1051 	         assert(varpos != varpos2);
1052 	
1053 	         /* Lapack uses only the lower left triangle of the symmetric matrix */
1054 	         quadmatrix[MIN(varpos, varpos2) * nexprs + MAX(varpos, varpos2)] = childcoefs[i] / 2.0;
1055 	      }
1056 	      else
1057 	      {
1058 	         varpos = SCIPhashmapGetImageInt(expr2idx, (void*) children[i]);
1059 	         assert(0 <= varpos && varpos < nexprs);
1060 	
1061 	         linvector[varpos] = childcoefs[i];
1062 	      }
1063 	   }
1064 	}
1065 	
1066 	/** tries to fill the nlhdlrexprdata for a potential quadratic SOC expression
1067 	 *
1068 	 * We say "try" because the expression might still turn out not to be a SOC at this point.
1069 	 */
1070 	static
1071 	SCIP_RETCODE tryFillNlhdlrExprDataQuad(
1072 	   SCIP*                 scip,               /**< SCIP data structure */
1073 	   SCIP_EXPR**           occurringexprs,     /**< array of all occurring expressions (nvars many) */
1074 	   SCIP_Real*            eigvecmatrix,       /**< array containing the Eigenvectors */
1075 	   SCIP_Real*            eigvals,            /**< array containing the Eigenvalues */
1076 	   SCIP_Real*            bp,                 /**< product of linear vector b * P (eigvecmatrix^t) */
1077 	   int                   nvars,              /**< number of variables */
1078 	   int*                  termbegins,         /**< pointer to store the termbegins */
1079 	   SCIP_Real*            transcoefs,         /**< pointer to store the transcoefs */
1080 	   int*                  transcoefsidx,      /**< pointer to store the transcoefsidx */
1081 	   SCIP_Real*            offsets,            /**< pointer to store the offsets */
1082 	   SCIP_Real*            lhsconstant,        /**< pointer to store the lhsconstant */
1083 	   int*                  nterms,             /**< pointer to store the total number of terms */
1084 	   SCIP_Bool*            success             /**< whether the expression is indeed a SOC */
1085 	   )
1086 	{
1087 	   SCIP_Real sqrteigval;
1088 	   int nextterm = 0;
1089 	   int nexttranscoef = 0;
1090 	   int specialtermidx;
1091 	   int i;
1092 	   int j;
1093 	
1094 	   assert(scip != NULL);
1095 	   assert(occurringexprs != NULL);
1096 	   assert(eigvecmatrix != NULL);
1097 	   assert(eigvals != NULL);
1098 	   assert(bp != NULL);
1099 	   assert(termbegins != NULL);
1100 	   assert(transcoefs != NULL);
1101 	   assert(transcoefsidx != NULL);
1102 	   assert(offsets != NULL);
1103 	   assert(lhsconstant != NULL);
1104 	   assert(success != NULL);
1105 	
1106 	   *success = FALSE;
1107 	   *nterms = 0;
1108 	
1109 	   /* we have lhsconstant + x^t A x + b x <= 0 and A has a single negative eigenvalue; try to build soc;
1110 	    * we now store all the v_i^T x + beta_i on the lhs, and compute the constant
1111 	    */
1112 	   specialtermidx = -1;
1113 	   for( i = 0; i < nvars; ++i )
1114 	   {
1115 	      if( SCIPisZero(scip, eigvals[i]) )
1116 	         continue;
1117 	
1118 	      if( eigvals[i] < 0.0 )
1119 	      {
1120 	         assert(specialtermidx == -1); /* there should only be one negative eigenvalue */
1121 	
1122 	         specialtermidx = i;
1123 	
1124 	         *lhsconstant -= bp[i] * bp[i] / (4.0 * eigvals[i]);
1125 	
1126 	         continue;
1127 	      }
1128 	
1129 	      assert(eigvals[i] > 0.0);
1130 	      sqrteigval = SQRT(eigvals[i]);
1131 	
1132 	      termbegins[nextterm] = nexttranscoef;
1133 	      offsets[nextterm] = bp[i] / (2.0 * sqrteigval);
1134 	      *lhsconstant -= bp[i] * bp[i] / (4.0 * eigvals[i]);
1135 	
1136 	      /* set transcoefs */
1137 	      for( j = 0; j < nvars; ++j )
1138 	      {
1139 	         if( !SCIPisZero(scip, eigvecmatrix[i * nvars + j]) )
1140 	         {
1141 	            transcoefs[nexttranscoef] = sqrteigval * eigvecmatrix[i * nvars + j];
1142 	            transcoefsidx[nexttranscoef] = j;
1143 	
1144 	            ++nexttranscoef;
1145 	         }
1146 	      }
1147 	      ++nextterm;
1148 	   }
1149 	   assert(specialtermidx > -1);
1150 	
1151 	   /* process constant; if constant is negative -> no soc */
1152 	   if( SCIPisNegative(scip, *lhsconstant) )
1153 	      return SCIP_OKAY;
1154 	
1155 	   /* we need lhsconstant to be >= 0 */
1156 	   if( *lhsconstant < 0.0 )
1157 	      *lhsconstant = 0.0;
1158 	
1159 	   /* store constant term */
1160 	   if( *lhsconstant > 0.0 )
1161 	   {
1162 	      termbegins[nextterm] = nexttranscoef;
1163 	      offsets[nextterm] = SQRT( *lhsconstant );
1164 	      ++nextterm;
1165 	   }
1166 	
1167 	   /* now process rhs */
1168 	   {
1169 	      SCIP_Real rhstermlb;
1170 	      SCIP_Real rhstermub;
1171 	      SCIP_Real signfactor;
1172 	
1173 	      assert(-eigvals[specialtermidx] > 0.0);
1174 	      sqrteigval = SQRT(-eigvals[specialtermidx]);
1175 	
1176 	      termbegins[nextterm] = nexttranscoef;
1177 	      offsets[nextterm] = -bp[specialtermidx] / (2.0 * sqrteigval);
1178 	
1179 	      /* the expression can only be an soc if the resulting rhs term does not change sign;
1180 	       * the rhs term is a linear combination of variables, so estimate its bounds
1181 	       */
1182 	      rhstermlb = offsets[nextterm];
1183 	      for( j = 0; j < nvars; ++j )
1184 	      {
1185 	         SCIP_INTERVAL activity;
1186 	         SCIP_Real aux;
1187 	
1188 	         if( SCIPisZero(scip, eigvecmatrix[specialtermidx * nvars + j]) )
1189 	            continue;
1190 	
1191 	         SCIP_CALL( SCIPevalExprActivity(scip, occurringexprs[j]) );
1192 	         activity = SCIPexprGetActivity(occurringexprs[j]);
1193 	
1194 	         if( eigvecmatrix[specialtermidx * nvars + j] > 0.0 )
1195 	         {
1196 	            aux = activity.inf;
1197 	            assert(!SCIPisInfinity(scip, aux));
1198 	         }
1199 	         else
1200 	         {
1201 	            aux = activity.sup;
1202 	            assert(!SCIPisInfinity(scip, -aux));
1203 	         }
1204 	
1205 	         if( SCIPisInfinity(scip, aux) || SCIPisInfinity(scip, -aux) )
1206 	         {
1207 	            rhstermlb = -SCIPinfinity(scip);
1208 	            break;
1209 	         }
1210 	         else
1211 	            rhstermlb += sqrteigval * eigvecmatrix[specialtermidx * nvars + j] * aux;
1212 	      }
1213 	
1214 	      rhstermub = offsets[nextterm];
1215 	      for( j = 0; j < nvars; ++j )
1216 	      {
1217 	         SCIP_INTERVAL activity;
1218 	         SCIP_Real aux;
1219 	
1220 	         if( SCIPisZero(scip, eigvecmatrix[specialtermidx * nvars + j]) )
1221 	            continue;
1222 	
1223 	         SCIP_CALL( SCIPevalExprActivity(scip, occurringexprs[j]) );
1224 	         activity = SCIPexprGetActivity(occurringexprs[j]);
1225 	
1226 	         if( eigvecmatrix[specialtermidx * nvars + j] > 0.0 )
1227 	         {
1228 	            aux = activity.sup;
1229 	            assert(!SCIPisInfinity(scip, -aux));
1230 	         }
1231 	         else
1232 	         {
1233 	            aux = activity.inf;
1234 	            assert(!SCIPisInfinity(scip, aux));
1235 	         }
1236 	
1237 	         if( SCIPisInfinity(scip, aux) || SCIPisInfinity(scip, -aux) )
1238 	         {
1239 	            rhstermub = SCIPinfinity(scip);
1240 	            break;
1241 	         }
1242 	         else
1243 	            rhstermub += sqrteigval * eigvecmatrix[specialtermidx * nvars + j] * aux;
1244 	      }
1245 	
1246 	      /* since we are just interested in obtaining an interval that contains the real bounds
1247 	       * and is tight enough so that we can identify that the rhsvar does not change sign,
1248 	       * we swap the bounds in case of numerical troubles
1249 	       */
1250 	      if( rhstermub < rhstermlb )
1251 	      {
1252 	         assert(SCIPisEQ(scip, rhstermub, rhstermlb));
1253 	         SCIPswapReals(&rhstermub, &rhstermlb);
1254 	      }
1255 	
1256 	      /* if rhs changes sign -> not a SOC */
1257 	      if( SCIPisLT(scip, rhstermlb, 0.0) && SCIPisGT(scip, rhstermub, 0.0) )
1258 	         return SCIP_OKAY;
1259 	
1260 	      signfactor = SCIPisLE(scip, rhstermub, 0.0) ? -1.0 : 1.0;
1261 	
1262 	      offsets[nextterm] *= signfactor;
1263 	
1264 	      /* set transcoefs for rhs term */
1265 	      for( j = 0; j < nvars; ++j )
1266 	      {
1267 	         if( SCIPisZero(scip, eigvecmatrix[specialtermidx * nvars + j]) )
1268 	            continue;
1269 	
1270 	         transcoefs[nexttranscoef] = signfactor * sqrteigval * eigvecmatrix[specialtermidx * nvars + j];
1271 	         transcoefsidx[nexttranscoef] = j;
1272 	
1273 	         ++nexttranscoef;
1274 	      }
1275 	
1276 	      /* if rhs is a constant this method shouldn't have been called */
1277 	      assert(nexttranscoef > termbegins[nextterm]);
1278 	
1279 	      /* finish processing term */
1280 	      ++nextterm;
1281 	   }
1282 	
1283 	   *nterms = nextterm;
1284 	
1285 	   /* sentinel value */
1286 	   termbegins[nextterm] = nexttranscoef;
1287 	
1288 	   *success = TRUE;
1289 	
1290 	   return SCIP_OKAY;
1291 	}
1292 	
1293 	/** detects if expr &le; auxvar is of the form SQRT(sum_i coef_i (expr_i + shift_i)^2 + const) &le; auxvar
1294 	 *
1295 	 * @note if a user inputs the above expression with `const` = -epsilon, then `const` is going to be set to 0.
1296 	 */
1297 	static
1298 	SCIP_RETCODE detectSocNorm(
1299 	   SCIP*                 scip,               /**< SCIP data structure */
1300 	   SCIP_EXPR*            expr,               /**< expression */
1301 	   SCIP_NLHDLREXPRDATA** nlhdlrexprdata,     /**< pointer to store nonlinear handler expression data */
1302 	   SCIP_Bool*            success             /**< pointer to store whether SOC structure has been detected */
1303 	   )
1304 	{
1305 	   SCIP_EXPR** children;
1306 	   SCIP_EXPR* child;
1307 	   SCIP_EXPR** vars;
1308 	   SCIP_HASHMAP* expr2idx;
1309 	   SCIP_HASHSET* linexprs;
1310 	   SCIP_Real* childcoefs;
1311 	   SCIP_Real* offsets;
1312 	   SCIP_Real* transcoefs;
1313 	   SCIP_Real constant;
1314 	   SCIP_Bool issoc;
1315 	   int* transcoefsidx;
1316 	   int* termbegins;
1317 	   int nchildren;
1318 	   int nterms;
1319 	   int nvars;
1320 	   int nextentry;
1321 	   int i;
1322 	
1323 	   assert(expr != NULL);
1324 	   assert(success != NULL);
1325 	
1326 	   *success = FALSE;
1327 	   issoc = TRUE;
1328 	
1329 	   /* relation is not "<=" -> skip */
1330 	   if( SCIPgetExprNLocksPosNonlinear(expr) == 0 )
1331 	      return SCIP_OKAY;
1332 	
1333 	   assert(SCIPexprGetNChildren(expr) > 0);
1334 	
1335 	   child = SCIPexprGetChildren(expr)[0];
1336 	   assert(child != NULL);
1337 	
1338 	   /* check whether expression is a SQRT and has a sum as child with at least 2 children and a non-negative constant */
1339 	   if( ! SCIPisExprPower(scip, expr)
1340 	      || SCIPgetExponentExprPow(expr) != 0.5
1341 	      || !SCIPisExprSum(scip, child)
1342 	      || SCIPexprGetNChildren(child) < 2
1343 	      || SCIPgetConstantExprSum(child) < 0.0)
1344 	   {
1345 	      return SCIP_OKAY;
1346 	   }
1347 	
1348 	   /* assert(SCIPvarGetLbLocal(auxvar) >= 0.0); */
1349 	
1350 	   /* get children of the sum */
1351 	   children = SCIPexprGetChildren(child);
1352 	   nchildren = SCIPexprGetNChildren(child);
1353 	   childcoefs = SCIPgetCoefsExprSum(child);
1354 	
1355 	   /* TODO: should we initialize the hashmap with size SCIPgetNVars() so that it never has to be resized? */
1356 	   SCIP_CALL( SCIPhashmapCreate(&expr2idx, SCIPblkmem(scip), nchildren) );
1357 	   SCIP_CALL( SCIPhashsetCreate(&linexprs, SCIPblkmem(scip), nchildren) );
1358 	
1359 	   /* we create coefs array here already, since we have to fill it in first loop in case of success
1360 	    * +1 for auxvar
1361 	    */
1362 	   SCIP_CALL( SCIPallocBufferArray(scip, &transcoefs, nchildren+1) );
1363 	
1364 	   nterms = 0;
1365 	
1366 	   /* check if all children are squares or linear terms with matching square term:
1367 	    * if the i-th child is (pow, expr, 2) we store the association <|expr -> i|> in expr2idx and if expr was in
1368 	    * linexprs, we remove it from there.
1369 	    * if the i-th child is expr' (different from (pow, expr, 2)) and expr' is not a key of expr2idx, we add it
1370 	    * to linexprs.
1371 	    * if at the end there is any expr in linexpr -> we do not have a separable quadratic function.
1372 	    */
1373 	   for( i = 0; i < nchildren; ++i )
1374 	   {
1375 	      /* handle quadratic expressions children */
1376 	      if( SCIPisExprPower(scip, children[i]) && SCIPgetExponentExprPow(children[i]) == 2.0 )
1377 	      {
1378 	         SCIP_EXPR* squarearg = SCIPexprGetChildren(children[i])[0];
1379 	
1380 	         if( !SCIPhashmapExists(expr2idx, (void*) squarearg) )
1381 	         {
1382 	            SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void *) squarearg, nterms) );
1383 	         }
1384 	
1385 	         if( childcoefs[i] < 0.0 )
1386 	         {
1387 	            issoc = FALSE;
1388 	            break;
1389 	         }
1390 	         transcoefs[nterms] = SQRT(childcoefs[i]);
1391 	
1392 	         SCIP_CALL( SCIPhashsetRemove(linexprs, (void*) squarearg) );
1393 	         ++nterms;
1394 	      }
1395 	      /* handle binary variable children */
1396 	      else if( SCIPisExprVar(scip, children[i]) && SCIPvarIsBinary(SCIPgetVarExprVar(children[i])) )
1397 	      {
1398 	         assert(!SCIPhashmapExists(expr2idx, (void*) children[i]));
1399 	         assert(!SCIPhashsetExists(linexprs, (void*) children[i]));
1400 	
1401 	         SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void *) children[i], nterms) );
1402 	
1403 	         if( childcoefs[i] < 0.0 )
1404 	         {
1405 	            issoc = FALSE;
1406 	            break;
1407 	         }
1408 	         transcoefs[nterms] = SQRT(childcoefs[i]);
1409 	
1410 	         ++nterms;
1411 	      }
1412 	      else
1413 	      {
1414 	         if( !SCIPhashmapExists(expr2idx, (void*) children[i]) )
1415 	         {
1416 	            SCIP_CALL( SCIPhashsetInsert(linexprs, SCIPblkmem(scip), (void*) children[i]) );
1417 	         }
1418 	      }
1419 	   }
1420 	
1421 	   /* there are linear terms without corresponding quadratic terms or it was detected not to be soc */
1422 	   if( SCIPhashsetGetNElements(linexprs) > 0 || ! issoc )
1423 	   {
1424 	      SCIPfreeBufferArray(scip, &transcoefs);
1425 	      SCIPhashsetFree(&linexprs, SCIPblkmem(scip) );
1426 	      SCIPhashmapFree(&expr2idx);
1427 	      return SCIP_OKAY;
1428 	   }
1429 	
1430 	   /* add one to terms counter for auxvar */
1431 	   ++nterms;
1432 	
1433 	   constant = SCIPgetConstantExprSum(child);
1434 	
1435 	   /* compute constant of possible soc expression to check its sign */
1436 	   for( i = 0; i < nchildren; ++i )
1437 	   {
1438 	      if( ! SCIPisExprPower(scip, children[i]) || SCIPgetExponentExprPow(children[i]) != 2.0 )
1439 	      {
1440 	         int auxvarpos;
1441 	
1442 	         assert(SCIPhashmapExists(expr2idx, (void*) children[i]) );
1443 	         auxvarpos = SCIPhashmapGetImageInt(expr2idx, (void*) children[i]);
1444 	
1445 	         constant -= SQR(0.5 * childcoefs[i] / transcoefs[auxvarpos]);
1446 	      }
1447 	   }
1448 	
1449 	   /* if the constant is negative -> no SOC */
1450 	   if( SCIPisNegative(scip, constant) )
1451 	   {
1452 	      SCIPfreeBufferArray(scip, &transcoefs);
1453 	      SCIPhashsetFree(&linexprs, SCIPblkmem(scip) );
1454 	      SCIPhashmapFree(&expr2idx);
1455 	      return SCIP_OKAY;
1456 	   }
1457 	   else if( SCIPisZero(scip, constant) )
1458 	      constant = 0.0;
1459 	   assert(constant >= 0.0);
1460 	
1461 	   /* at this point, we have found an SOC structure */
1462 	   *success = TRUE;
1463 	
1464 	   nvars = nterms;
1465 	
1466 	   /* add one to terms counter for constant term */
1467 	   if( constant > 0.0 )
1468 	      ++nterms;
1469 	
1470 	   /* allocate temporary memory to collect data */
1471 	   SCIP_CALL( SCIPallocBufferArray(scip, &vars, nvars) );
1472 	   SCIP_CALL( SCIPallocBufferArray(scip, &offsets, nterms) );
1473 	   SCIP_CALL( SCIPallocBufferArray(scip, &transcoefsidx, nvars) );
1474 	   SCIP_CALL( SCIPallocBufferArray(scip, &termbegins, nterms + 1) );
1475 	
1476 	   /* fill in data for non constant terms of lhs; initialize their offsets */
1477 	   for( i = 0; i < nvars - 1; ++i )
1478 	   {
1479 	      transcoefsidx[i] = i;
1480 	      termbegins[i] = i;
1481 	      offsets[i] = 0.0;
1482 	   }
1483 	
1484 	   /* add constant term and rhs */
1485 	   vars[nvars - 1] = expr;
1486 	   if( constant > 0.0 )
1487 	   {
1488 	      /* constant term */
1489 	      termbegins[nterms - 2] = nterms - 2;
1490 	      offsets[nterms - 2] = SQRT(constant);
1491 	
1492 	      /* rhs */
1493 	      termbegins[nterms - 1] = nterms - 2;
1494 	      offsets[nterms - 1] = 0.0;
1495 	      transcoefsidx[nterms - 2] = nvars - 1;
1496 	      transcoefs[nterms - 2] = 1.0;
1497 	
1498 	      /* sentinel value */
1499 	      termbegins[nterms] = nterms - 1;
1500 	   }
1501 	   else
1502 	   {
1503 	      /* rhs */
1504 	      termbegins[nterms - 1] = nterms - 1;
1505 	      offsets[nterms - 1] = 0.0;
1506 	      transcoefsidx[nterms - 1] = nvars - 1;
1507 	      transcoefs[nterms - 1] = 1.0;
1508 	
1509 	      /* sentinel value */
1510 	      termbegins[nterms] = nterms;
1511 	   }
1512 	
1513 	   /* request required auxiliary variables and fill vars and offsets array */
1514 	   nextentry = 0;
1515 	   for( i = 0; i < nchildren; ++i )
1516 	   {
1517 	      if( SCIPisExprPower(scip, children[i]) && SCIPgetExponentExprPow(children[i]) == 2.0 )
1518 	      {
1519 	         SCIP_EXPR* squarearg;
1520 	
1521 	         squarearg = SCIPexprGetChildren(children[i])[0];
1522 	         assert(SCIPhashmapGetImageInt(expr2idx, (void*) squarearg) == nextentry);
1523 	
1524 	         SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, squarearg, TRUE, FALSE, FALSE, FALSE) );
1525 	
1526 	         vars[nextentry] = squarearg;
1527 	         ++nextentry;
1528 	      }
1529 	      else if( SCIPisExprVar(scip, children[i]) && SCIPvarIsBinary(SCIPgetVarExprVar(children[i])) )
1530 	      {
1531 	         /* handle binary variable children: no need to request auxvar */
1532 	         assert(SCIPhashmapGetImageInt(expr2idx, (void*) children[i]) == nextentry);
1533 	         vars[nextentry] = children[i];
1534 	         ++nextentry;
1535 	      }
1536 	      else
1537 	      {
1538 	         int auxvarpos;
1539 	
1540 	         assert(SCIPhashmapExists(expr2idx, (void*) children[i]));
1541 	         auxvarpos = SCIPhashmapGetImageInt(expr2idx, (void*) children[i]);
1542 	
1543 	         SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, children[i], TRUE, FALSE, FALSE, FALSE) );
1544 	
1545 	         offsets[auxvarpos] = 0.5 * childcoefs[i] / transcoefs[auxvarpos];
1546 	      }
1547 	   }
1548 	   assert(nextentry == nvars - 1);
1549 	
1550 	#ifdef SCIP_DEBUG
1551 	   SCIPdebugMsg(scip, "found SOC structure for expression %p\n", (void*)expr);
1552 	   SCIPprintExpr(scip, expr, NULL);
1553 	   SCIPinfoMessage(scip, NULL, " <= auxvar\n");
1554 	#endif
1555 	
1556 	   /* create and store nonlinear handler expression data */
1557 	   SCIP_CALL( createNlhdlrExprData(scip, vars, offsets, transcoefs, transcoefsidx, termbegins,
1558 	            nvars, nterms, nlhdlrexprdata) );
1559 	   assert(*nlhdlrexprdata != NULL);
1560 	
1561 	   /* free memory */
1562 	   SCIPhashsetFree(&linexprs, SCIPblkmem(scip) );
1563 	   SCIPhashmapFree(&expr2idx);
1564 	   SCIPfreeBufferArray(scip, &termbegins);
1565 	   SCIPfreeBufferArray(scip, &transcoefsidx);
1566 	   SCIPfreeBufferArray(scip, &offsets);
1567 	   SCIPfreeBufferArray(scip, &vars);
1568 	   SCIPfreeBufferArray(scip, &transcoefs);
1569 	
1570 	   return SCIP_OKAY;
1571 	}
1572 	
1573 	/** helper method to detect c + sum_i coef_i expr_i^2 - coef_k expr_k^2 &le; 0
1574 	 *  and c + sum_i coef_i expr_i^2 - coef_k expr_k expr_l &le; 0
1575 	 *
1576 	 *  binary linear variables are interpreted as quadratic terms
1577 	 *
1578 	 *  @todo: extend this function to detect  c + sum_i coef_i (expr_i + const_i)^2 - ...
1579 	 *  this would probably share a lot of code with detectSocNorm
1580 	 */
1581 	static
1582 	SCIP_RETCODE detectSocQuadraticSimple(
1583 	   SCIP*                 scip,               /**< SCIP data structure */
1584 	   SCIP_EXPR*            expr,               /**< expression */
1585 	   SCIP_Real             conslhs,            /**< lhs of the constraint that the expression defines (or SCIP_INVALID) */
1586 	   SCIP_Real             consrhs,            /**< rhs of the constraint that the expression defines (or SCIP_INVALID) */
1587 	   SCIP_NLHDLREXPRDATA** nlhdlrexprdata,     /**< pointer to store nonlinear handler expression data */
1588 	   SCIP_Bool*            enforcebelow,       /**< pointer to store whether we enforce <= (TRUE) or >= (FALSE); only valid when success is TRUE */
1589 	   SCIP_Bool*            success             /**< pointer to store whether SOC structure has been detected */
1590 	   )
1591 	{
1592 	   SCIP_EXPR** children;
1593 	   SCIP_EXPR** vars = NULL;
1594 	   SCIP_Real* childcoefs;
1595 	   SCIP_Real* offsets = NULL;
1596 	   SCIP_Real* transcoefs = NULL;
1597 	   int* transcoefsidx = NULL;
1598 	   int* termbegins = NULL;
1599 	   SCIP_Real constant;
1600 	   SCIP_Real lhsconstant;
1601 	   SCIP_Real lhs;
1602 	   SCIP_Real rhs;
1603 	   SCIP_Real rhssign;
1604 	   SCIP_INTERVAL expractivity;
1605 	   int ntranscoefs;
1606 	   int nposquadterms;
1607 	   int nnegquadterms;
1608 	   int nposbilinterms;
1609 	   int nnegbilinterms;
1610 	   int rhsidx;
1611 	   int lhsidx;
1612 	   int specialtermidx;
1613 	   int nchildren;
1614 	   int nnzinterms;
1615 	   int nterms;
1616 	   int nvars;
1617 	   int nextentry;
1618 	   int i;
1619 	   SCIP_Bool ishyperbolic;
1620 	
1621 	   assert(expr != NULL);
1622 	   assert(success != NULL);
1623 	
1624 	   *success = FALSE;
1625 	
1626 	   /* check whether expression is a sum with at least 2 children */
1627 	   if( ! SCIPisExprSum(scip, expr) || SCIPexprGetNChildren(expr) < 2 )
1628 	      return SCIP_OKAY;
1629 	
1630 	   /* get children of the sum */
1631 	   children = SCIPexprGetChildren(expr);
1632 	   nchildren = SCIPexprGetNChildren(expr);
1633 	   constant = SCIPgetConstantExprSum(expr);
1634 	
1635 	   /* we duplicate the child coefficients since we have to manipulate them */
1636 	   SCIP_CALL( SCIPduplicateBufferArray(scip, &childcoefs, SCIPgetCoefsExprSum(expr), nchildren) ); /*lint !e666*/
1637 	
1638 	   /* initialize data */
1639 	   lhsidx = -1;
1640 	   rhsidx = -1;
1641 	   nposquadterms = 0;
1642 	   nnegquadterms = 0;
1643 	   nposbilinterms = 0;
1644 	   nnegbilinterms = 0;
1645 	
1646 	   /* check if all children are quadratic or binary linear and count number of positive and negative terms */
1647 	   for( i = 0; i < nchildren; ++i )
1648 	   {
1649 	      if( SCIPisExprPower(scip, children[i]) && SCIPgetExponentExprPow(children[i]) == 2.0 )
1650 	      {
1651 	         if( childcoefs[i] > 0.0 )
1652 	         {
1653 	            ++nposquadterms;
1654 	            lhsidx = i;
1655 	         }
1656 	         else
1657 	         {
1658 	            ++nnegquadterms;
1659 	            rhsidx = i;
1660 	         }
1661 	      }
1662 	      else if( SCIPisExprVar(scip, children[i]) && SCIPvarIsBinary(SCIPgetVarExprVar(children[i])) )
1663 	      {
1664 	         if( childcoefs[i] > 0.0 )
1665 	         {
1666 	            ++nposquadterms;
1667 	            lhsidx = i;
1668 	         }
1669 	         else
1670 	         {
1671 	            ++nnegquadterms;
1672 	            rhsidx = i;
1673 	         }
1674 	      }
1675 	      else if( SCIPisExprProduct(scip, children[i]) && SCIPexprGetNChildren(children[i]) == 2 )
1676 	      {
1677 	         if( childcoefs[i] > 0.0 )
1678 	         {
1679 	            ++nposbilinterms;
1680 	            lhsidx = i;
1681 	         }
1682 	         else
1683 	         {
1684 	            ++nnegbilinterms;
1685 	            rhsidx = i;
1686 	         }
1687 	      }
1688 	      else
1689 	      {
1690 	         goto CLEANUP;
1691 	      }
1692 	
1693 	      /* more than one positive eigenvalue and more than one negative eigenvalue -> can't be convex */
1694 	      if( nposquadterms > 1 && nnegquadterms > 1 )
1695 	         goto CLEANUP;
1696 	
1697 	      /* more than one bilinear term -> can't be handled by this method */
1698 	      if( nposbilinterms + nnegbilinterms > 1 )
1699 	         goto CLEANUP;
1700 	
1701 	      /* one positive bilinear term and also at least one positive quadratic term -> not a simple SOC */
1702 	      if( nposbilinterms > 0 && nposquadterms > 0 )
1703 	         goto CLEANUP;
1704 	
1705 	      /* one negative bilinear term and also at least one negative quadratic term -> not a simple SOC */
1706 	      if( nnegbilinterms > 0 && nnegquadterms > 0 )
1707 	         goto CLEANUP;
1708 	   }
1709 	
1710 	   if( nposquadterms == nchildren || nnegquadterms == nchildren )
1711 	      goto CLEANUP;
1712 	
1713 	   assert(nposquadterms <= 1 || nnegquadterms <= 1);
1714 	   assert(nposbilinterms + nnegbilinterms <= 1);
1715 	   assert(nposbilinterms == 0 || nposquadterms == 0);
1716 	   assert(nnegbilinterms == 0 || nnegquadterms == 0);
1717 	
1718 	   /* if a bilinear term is involved, it is a hyperbolic expression */
1719 	   ishyperbolic = (nposbilinterms + nnegbilinterms > 0);
1720 	
1721 	   if( conslhs == SCIP_INVALID || consrhs == SCIP_INVALID )  /*lint !e777*/
1722 	   {
1723 	      SCIP_CALL( SCIPevalExprActivity(scip, expr) );
1724 	      expractivity = SCIPexprGetActivity(expr);
1725 	
1726 	      lhs = (conslhs == SCIP_INVALID ? expractivity.inf : conslhs); /*lint !e777*/
1727 	      rhs = (consrhs == SCIP_INVALID ? expractivity.sup : consrhs); /*lint !e777*/
1728 	   }
1729 	   else
1730 	   {
1731 	      lhs = conslhs;
1732 	      rhs = consrhs;
1733 	   }
1734 	
1735 	   /* detect case and store lhs/rhs information */
1736 	   if( (ishyperbolic && nnegbilinterms > 0) || (!ishyperbolic && nnegquadterms < 2) )
1737 	   {
1738 	      /* we have -x*y + z^2 ... -> we want to write  z^2 ... <= x*y;
1739 	       * or we have -x^2 + y^2  ... -> we want to write y^2 ... <= x^2;
1740 	       * in any case, we need a finite rhs
1741 	       */
1742 	      assert(nnegbilinterms == 1 || nnegquadterms == 1);
1743 	      assert(rhsidx != -1);
1744 	
1745 	      /* if rhs is infinity, it can't be soc
1746 	       * TODO: if it can't be soc, then we should enforce the caller so that we do not try the more complex quadratic
1747 	       * method
1748 	       */
1749 	      if( SCIPisInfinity(scip, rhs) )
1750 	         goto CLEANUP;
1751 	
1752 	      specialtermidx = rhsidx;
1753 	      lhsconstant = constant - rhs;
1754 	      *enforcebelow = TRUE; /* enforce expr <= rhs */
1755 	   }
1756 	   else
1757 	   {
1758 	      assert(lhsidx != -1);
1759 	
1760 	      /* if lhs is infinity, it can't be soc */
1761 	      if( SCIPisInfinity(scip, -lhs) )
1762 	         goto CLEANUP;
1763 	
1764 	      specialtermidx = lhsidx;
1765 	      lhsconstant = lhs - constant;
1766 	
1767 	      /* negate all coefficients */
1768 	      for( i = 0; i < nchildren; ++i )
1769 	         childcoefs[i] = -childcoefs[i];
1770 	      *enforcebelow = FALSE; /* enforce lhs <= expr */
1771 	   }
1772 	   assert(childcoefs[specialtermidx] != 0.0);
1773 	
1774 	   if( ishyperbolic )
1775 	   {
1776 	      SCIP_INTERVAL yactivity;
1777 	      SCIP_INTERVAL zactivity;
1778 	
1779 	      assert(SCIPexprGetNChildren(children[specialtermidx]) == 2);
1780 	
1781 	      SCIP_CALL( SCIPevalExprActivity(scip, SCIPexprGetChildren(children[specialtermidx])[0]) );
1782 	      yactivity = SCIPexprGetActivity(SCIPexprGetChildren(children[specialtermidx])[0]);
1783 	
1784 	      SCIP_CALL( SCIPevalExprActivity(scip, SCIPexprGetChildren(children[specialtermidx])[1]) );
1785 	      zactivity = SCIPexprGetActivity(SCIPexprGetChildren(children[specialtermidx])[1]);
1786 	
1787 	      if( SCIPisNegative(scip, yactivity.inf + zactivity.inf) )
1788 	      {
1789 	         /* the sum of the expressions in the bilinear term changes sign -> no SOC */
1790 	         if( SCIPisPositive(scip, yactivity.sup + zactivity.sup) )
1791 	            goto CLEANUP;
1792 	
1793 	         rhssign = -1.0;
1794 	      }
1795 	      else
1796 	         rhssign = 1.0;
1797 	
1798 	      lhsconstant *= 4.0 / -childcoefs[specialtermidx];
1799 	   }
1800 	   else if( SCIPisExprVar(scip, children[specialtermidx]) )
1801 	   {
1802 	      /* children[specialtermidx] can be a variable, in which case we treat it as if it is squared */
1803 	      rhssign = 1.0;
1804 	   }
1805 	   else
1806 	   {
1807 	      SCIP_INTERVAL rhsactivity;
1808 	
1809 	      assert(SCIPexprGetNChildren(children[specialtermidx]) == 1);
1810 	      SCIP_CALL( SCIPevalExprActivity(scip, SCIPexprGetChildren(children[specialtermidx])[0]) );
1811 	      rhsactivity = SCIPexprGetActivity(SCIPexprGetChildren(children[specialtermidx])[0]);
1812 	
1813 	      if( rhsactivity.inf < 0.0 )
1814 	      {
1815 	         /* rhs variable changes sign -> no SOC */
1816 	         if( rhsactivity.sup > 0.0 )
1817 	            goto CLEANUP;
1818 	
1819 	         rhssign = -1.0;
1820 	      }
1821 	      else
1822 	         rhssign = 1.0;
1823 	   }
1824 	
1825 	   if( SCIPisNegative(scip, lhsconstant) )
1826 	      goto CLEANUP;
1827 	
1828 	   if( SCIPisZero(scip, lhsconstant) )
1829 	      lhsconstant = 0.0;
1830 	
1831 	   /*
1832 	    * we have found an SOC-representable expression. Now build the nlhdlrexprdata
1833 	    *
1834 	    * in the non-hyperbolic case, c + sum_i coef_i expr_i^2 - coef_k expr_k^2 <= 0 is transformed to
1835 	    * SQRT( c + sum_i coef_i expr_i^2 ) <= coef_k expr_k
1836 	    * so there are nchildren many vars, nchildren (+ 1 if c != 0) many terms, nchildren many coefficients in the vs
1837 	    * in SOC representation
1838 	    *
1839 	    * in the hyperbolic case, c + sum_i coef_i expr_i^2 - coef_k expr_k expr_l <= 0 is transformed to
1840 	    * SQRT( 4(c + sum_i coef_i expr_i^2) + (expr_k - expr_l)^2 ) <= expr_k + expr_l
1841 	    * so there are nchildren + 1many vars, nchildren + 1(+ 1 if c != 0) many terms, nchildren + 3 many coefficients in
1842 	    * the vs in SOC representation
1843 	    */
1844 	
1845 	   ntranscoefs = ishyperbolic ? nchildren + 3 : nchildren;
1846 	   nvars = ishyperbolic ? nchildren + 1 : nchildren;
1847 	   nterms = nvars;
1848 	
1849 	   /* constant term */
1850 	   if( lhsconstant > 0.0 )
1851 	      nterms++;
1852 	
1853 	   /* SOC was detected, allocate temporary memory for data to collect */
1854 	   SCIP_CALL( SCIPallocBufferArray(scip, &vars, nvars) );
1855 	   SCIP_CALL( SCIPallocBufferArray(scip, &offsets, nterms) );
1856 	   SCIP_CALL( SCIPallocBufferArray(scip, &transcoefs, ntranscoefs) );
1857 	   SCIP_CALL( SCIPallocBufferArray(scip, &transcoefsidx, ntranscoefs) );
1858 	   SCIP_CALL( SCIPallocBufferArray(scip, &termbegins, nterms + 1) );
1859 	
1860 	   *success = TRUE;
1861 	   nextentry = 0;
1862 	
1863 	   /* collect all the v_i and beta_i */
1864 	   nnzinterms = 0;
1865 	   for( i = 0; i < nchildren; ++i )
1866 	   {
1867 	      /* variable and coef for rhs have to be set to the last entry */
1868 	      if( i == specialtermidx )
1869 	         continue;
1870 	
1871 	      /* extract (unique) variable appearing in term */
1872 	      if( SCIPisExprVar(scip, children[i]) )
1873 	      {
1874 	         vars[nextentry] = children[i];
1875 	
1876 	         assert(SCIPvarIsBinary(SCIPgetVarExprVar(vars[nextentry])));
1877 	      }
1878 	      else
1879 	      {
1880 	         assert(SCIPisExprPower(scip, children[i]));
1881 	
1882 	         /* notify that we will require auxiliary variable */
1883 	         SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, SCIPexprGetChildren(children[i])[0], TRUE, FALSE, FALSE, FALSE) );
1884 	         vars[nextentry] = SCIPexprGetChildren(children[i])[0];
1885 	      }
1886 	      assert(vars[nextentry] != NULL);
1887 	
1888 	      /* store v_i and beta_i */
1889 	      termbegins[nextentry] = nnzinterms;
1890 	      offsets[nextentry] = 0.0;
1891 	
1892 	      transcoefsidx[nnzinterms] = nextentry;
1893 	      if( ishyperbolic )
1894 	      {
1895 	         /* we eliminate the coefficient of the bilinear term to arrive at standard form */
1896 	         assert(4.0 * childcoefs[i] / -childcoefs[specialtermidx] > 0.0);
1897 	         transcoefs[nnzinterms] = SQRT(4.0 * childcoefs[i] / -childcoefs[specialtermidx]);
1898 	      }
1899 	      else
1900 	      {
1901 	         assert(childcoefs[i] > 0.0);
1902 	         transcoefs[nnzinterms] = SQRT(childcoefs[i]);
1903 	      }
1904 	
1905 	      /* finish adding nonzeros */
1906 	      ++nnzinterms;
1907 	
1908 	      /* finish processing term */
1909 	      ++nextentry;
1910 	   }
1911 	   assert(nextentry == nchildren - 1);
1912 	
1913 	   /* store term for constant (v_i = 0) */
1914 	   if( lhsconstant > 0.0 )
1915 	   {
1916 	      termbegins[nextentry] = nnzinterms;
1917 	      offsets[nextentry] = SQRT(lhsconstant);
1918 	
1919 	      /* finish processing term; this term has 0 nonzero thus we do not increase nnzinterms */
1920 	      ++nextentry;
1921 	   }
1922 	
1923 	   if( !ishyperbolic )
1924 	   {
1925 	      /* store rhs term */
1926 	      if( SCIPisExprVar(scip, children[specialtermidx]) )
1927 	      {
1928 	         /* this should be the "children[specialtermidx] can be a variable, in which case we treat it as if it is squared" case */
1929 	         SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, children[specialtermidx], TRUE, FALSE, FALSE, FALSE) );
1930 	         vars[nvars - 1] = children[specialtermidx];
1931 	      }
1932 	      else
1933 	      {
1934 	         assert(SCIPisExprPower(scip, children[specialtermidx]));
1935 	         assert(SCIPexprGetChildren(children[specialtermidx]) != NULL);
1936 	         SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, SCIPexprGetChildren(children[specialtermidx])[0], TRUE, FALSE, FALSE, FALSE) );
1937 	         vars[nvars - 1] = SCIPexprGetChildren(children[specialtermidx])[0];
1938 	      }
1939 	
1940 	      assert(childcoefs[specialtermidx] < 0.0);
1941 	
1942 	      termbegins[nextentry] = nnzinterms;
1943 	      offsets[nextentry] = 0.0;
1944 	      transcoefs[nnzinterms] = rhssign * SQRT(-childcoefs[specialtermidx]);
1945 	      transcoefsidx[nnzinterms] = nvars - 1;
1946 	
1947 	      /* finish adding nonzeros */
1948 	      ++nnzinterms;
1949 	
1950 	      /* finish processing term */
1951 	      ++nextentry;
1952 	   }
1953 	   else
1954 	   {
1955 	      /* store last lhs term and rhs term coming from the bilinear term */
1956 	      SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, SCIPexprGetChildren(children[specialtermidx])[0], TRUE, FALSE, FALSE, FALSE) );
1957 	      vars[nvars - 2] = SCIPexprGetChildren(children[specialtermidx])[0];
1958 	
1959 	      SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, SCIPexprGetChildren(children[specialtermidx])[1], TRUE, FALSE, FALSE, FALSE) );
1960 	      vars[nvars - 1] = SCIPexprGetChildren(children[specialtermidx])[1];
1961 	
1962 	      /* at this point, vars[nvars - 2] = expr_k and vars[nvars - 1] = expr_l;
1963 	       * on the lhs we have the term (expr_k - expr_l)^2
1964 	       */
1965 	      termbegins[nextentry] = nnzinterms;
1966 	      offsets[nextentry] = 0.0;
1967 	
1968 	      /* expr_k */
1969 	      transcoefsidx[nnzinterms] = nvars - 2;
1970 	      transcoefs[nnzinterms] = 1.0;
1971 	      ++nnzinterms;
1972 	
1973 	      /* - expr_l */
1974 	      transcoefsidx[nnzinterms] = nvars - 1;
1975 	      transcoefs[nnzinterms] = -1.0;
1976 	      ++nnzinterms;
1977 	
1978 	      /* finish processing term */
1979 	      ++nextentry;
1980 	
1981 	      /* on rhs we have +/-(expr_k + expr_l) */
1982 	      termbegins[nextentry] = nnzinterms;
1983 	      offsets[nextentry] = 0.0;
1984 	
1985 	      /* rhssing * expr_k */
1986 	      transcoefsidx[nnzinterms] = nvars - 2;
1987 	      transcoefs[nnzinterms] = rhssign;
1988 	      ++nnzinterms;
1989 	
1990 	      /* rhssing * expr_l */
1991 	      transcoefsidx[nnzinterms] = nvars - 1;
1992 	      transcoefs[nnzinterms] = rhssign;
1993 	      ++nnzinterms;
1994 	
1995 	      /* finish processing term */
1996 	      ++nextentry;
1997 	   }
1998 	   assert(nextentry == nterms);
1999 	   assert(nnzinterms == ntranscoefs);
2000 	
2001 	   /* sentinel value */
2002 	   termbegins[nextentry] = nnzinterms;
2003 	
2004 	#ifdef SCIP_DEBUG
2005 	   SCIPdebugMsg(scip, "found SOC structure for expression %p\n  %g <= ", (void*)expr, lhs);
2006 	   SCIPprintExpr(scip, expr, NULL);
2007 	   SCIPinfoMessage(scip, NULL, " <= %g\n", rhs);
2008 	#endif
2009 	
2010 	   /* create and store nonlinear handler expression data */
2011 	   SCIP_CALL( createNlhdlrExprData(scip, vars, offsets, transcoefs, transcoefsidx, termbegins, nvars, nterms,
2012 	            nlhdlrexprdata) );
2013 	   assert(*nlhdlrexprdata != NULL);
2014 	
2015 	CLEANUP:
2016 	   SCIPfreeBufferArrayNull(scip, &termbegins);
2017 	   SCIPfreeBufferArrayNull(scip, &transcoefsidx);
2018 	   SCIPfreeBufferArrayNull(scip, &transcoefs);
2019 	   SCIPfreeBufferArrayNull(scip, &offsets);
2020 	   SCIPfreeBufferArrayNull(scip, &vars);
2021 	   SCIPfreeBufferArrayNull(scip, &childcoefs);
2022 	
2023 	   return SCIP_OKAY;
2024 	}
2025 	
2026 	/** detects complex quadratic expressions that can be represented as SOC constraints
2027 	 *
2028 	 *  These are quadratic expressions with either exactly one positive or exactly one negative eigenvalue,
2029 	 *  in addition to some extra conditions. One needs to write the quadratic as
2030 	 *  sum eigval_i (eigvec_i . x)^2 + c &le; -eigval_k (eigvec_k . x)^2, where eigval_k is the negative eigenvalue,
2031 	 *  and c must be positive and (eigvec_k . x) must not change sign.
2032 	 *  This is described in more details in
2033 	 *  Mahajan, Ashutosh & Munson, Todd, Exploiting Second-Order Cone Structure for Global Optimization, 2010.
2034 	 *
2035 	 *  The eigen-decomposition is computed using Lapack.
2036 	 *  Binary linear variables are interpreted as quadratic terms.
2037 	 *
2038 	 * @todo: In the case -b <= a + x^2 - y^2 <= b, it is possible to represent both sides by SOC. Currently, the
2039 	 * datastructure can only handle one SOC. If this should appear more often, it could be worth to extend it,
2040 	 * such that both sides can be handled (see e.g. instance chp_partload).
2041 	 * FS: this shouldn't be possible. For a <= b + x^2 - y^2 <= c to be SOC representable on both sides, we would need
2042 	 * that a - b >= 0 and b -c >= 0, but this implies that a >= c and assuming the constraint is not trivially infeasible,
2043 	 * a <= b. Thus, a = b = c and the constraint is x^2 == y^2.
2044 	 *
2045 	 * @todo: Since cons_nonlinear multiplies as many terms out as possible during presolving, some SOC-representable
2046 	 * structures cannot be detected, (see e.g. instances bearing or wager). There is currently no obvious way
2047 	 * to handle this.
2048 	 */
2049 	static
2050 	SCIP_RETCODE detectSocQuadraticComplex(
2051 	   SCIP*                 scip,               /**< SCIP data structure */
2052 	   SCIP_EXPR*            expr,               /**< expression */
2053 	   SCIP_Real             conslhs,            /**< lhs of the constraint that the expression defines (or SCIP_INVALID) */
2054 	   SCIP_Real             consrhs,            /**< rhs of the constraint that the expression defines (or SCIP_INVALID) */
2055 	   SCIP_NLHDLREXPRDATA** nlhdlrexprdata,     /**< pointer to store nonlinear handler expression data */
2056 	   SCIP_Bool*            enforcebelow,       /**< pointer to store whether we enforce <= (TRUE) or >= (FALSE); only
2057 	                                              *   valid when success is TRUE */
2058 	   SCIP_Bool*            success             /**< pointer to store whether SOC structure has been detected */
2059 	   )
2060 	{
2061 	   SCIP_EXPR** occurringexprs;
2062 	   SCIP_HASHMAP* expr2idx;
2063 	   SCIP_Real* offsets;
2064 	   SCIP_Real* transcoefs;
2065 	   SCIP_Real* eigvecmatrix;
2066 	   SCIP_Real* eigvals;
2067 	   SCIP_Real* lincoefs;
2068 	   SCIP_Real* bp;
2069 	   int* transcoefsidx;
2070 	   int* termbegins;
2071 	   SCIP_Real constant;
2072 	   SCIP_Real lhsconstant;
2073 	   SCIP_Real lhs;
2074 	   SCIP_Real rhs;
2075 	   SCIP_INTERVAL expractivity;
2076 	   int nvars;
2077 	   int nterms;
2078 	   int nchildren;
2079 	   int npos;
2080 	   int nneg;
2081 	   int ntranscoefs;
2082 	   int i;
2083 	   int j;
2084 	   SCIP_Bool rhsissoc;
2085 	   SCIP_Bool lhsissoc;
2086 	   SCIP_Bool isquadratic;
2087 	
2088 	   assert(expr != NULL);
2089 	   assert(success != NULL);
2090 	
2091 	   *success = FALSE;
2092 	
2093 	   /* check whether expression is a sum with at least 2 children */
2094 	   if( ! SCIPisExprSum(scip, expr) || SCIPexprGetNChildren(expr) < 2 )
2095 	   {
2096 	      return SCIP_OKAY;
2097 	   }
2098 	
2099 	   /* we need Lapack to compute eigenvalues/vectors below */
2100 	   if( ! SCIPlapackIsAvailable() )
2101 	      return SCIP_OKAY;
2102 	
2103 	   /* get children of the sum */
2104 	   nchildren = SCIPexprGetNChildren(expr);
2105 	   constant = SCIPgetConstantExprSum(expr);
2106 	
2107 	   /* initialize data */
2108 	   offsets = NULL;
2109 	   transcoefs = NULL;
2110 	   transcoefsidx = NULL;
2111 	   termbegins = NULL;
2112 	   bp = NULL;
2113 	
2114 	   SCIP_CALL( SCIPhashmapCreate(&expr2idx, SCIPblkmem(scip), 2 * nchildren) );
2115 	   SCIP_CALL( SCIPallocBufferArray(scip, &occurringexprs, 2 * nchildren) );
2116 	
2117 	   /* check if the expression is quadratic and collect all occurring expressions */
2118 	   SCIP_CALL( checkAndCollectQuadratic(scip, expr, expr2idx, occurringexprs, &nvars, &isquadratic) );
2119 	
2120 	   if( !isquadratic )
2121 	   {
2122 	      SCIPfreeBufferArray(scip, &occurringexprs);
2123 	      SCIPhashmapFree(&expr2idx);
2124 	      return SCIP_OKAY;
2125 	   }
2126 	
2127 	   /* check that nvars*nvars doesn't get too large, see also SCIPcomputeExprQuadraticCurvature() */
2128 	   if( nvars > 7000 )
2129 	   {
2130 	      SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "nlhdlr_soc - number of quadratic variables is too large (%d) to check the curvature\n", nvars);
2131 	      SCIPfreeBufferArray(scip, &occurringexprs);
2132 	      SCIPhashmapFree(&expr2idx);
2133 	      return SCIP_OKAY;
2134 	   }
2135 	
2136 	   assert(SCIPhashmapGetNElements(expr2idx) == nvars);
2137 	
2138 	   /* create datastructures for constaint defining matrix and vector */
2139 	   SCIP_CALL( SCIPallocClearBufferArray(scip, &eigvecmatrix, nvars * nvars) ); /*lint !e647*/
2140 	   SCIP_CALL( SCIPallocClearBufferArray(scip, &lincoefs, nvars) );
2141 	
2142 	   /* build constraint defining matrix (stored in eigvecmatrix) and vector (stored in lincoefs) */
2143 	   buildQuadExprMatrix(scip, expr, expr2idx, nvars, eigvecmatrix, lincoefs);
2144 	
2145 	   SCIP_CALL( SCIPallocBufferArray(scip, &eigvals, nvars) );
2146 	
2147 	   /* compute eigenvalues and vectors, A = PDP^t
2148 	    * note: eigvecmatrix stores P^t, i.e., P^t_{i,j} = eigvecmatrix[i*nvars+j]
2149 	    */
2150 	   if( SCIPlapackComputeEigenvalues(SCIPbuffer(scip), TRUE, nvars, eigvecmatrix, eigvals) != SCIP_OKAY )
2151 	   {
2152 	      SCIPdebugMsg(scip, "Failed to compute eigenvalues and eigenvectors for expression:\n");
2153 	
2154 	#ifdef SCIP_DEBUG
2155 	      SCIPdismantleExpr(scip, NULL, expr);
2156 	#endif
2157 	
2158 	      goto CLEANUP;
2159 	   }
2160 	
2161 	   SCIP_CALL( SCIPallocClearBufferArray(scip, &bp, nvars) );
2162 	
2163 	   nneg = 0;
2164 	   npos = 0;
2165 	   ntranscoefs = 0;
2166 	
2167 	   /* set small eigenvalues to 0 and compute b*P */
2168 	   for( i = 0; i < nvars; ++i )
2169 	   {
2170 	      for( j = 0; j < nvars; ++j )
2171 	      {
2172 	         bp[i] += lincoefs[j] * eigvecmatrix[i * nvars + j];
2173 	
2174 	         /* count the number of transcoefs to be used later */
2175 	         if( !SCIPisZero(scip, eigvals[i]) && !SCIPisZero(scip, eigvecmatrix[i * nvars + j]) )
2176 	            ++ntranscoefs;
2177 	      }
2178 	
2179 	      if( SCIPisZero(scip, eigvals[i]) )
2180 	      {
2181 	         /* if there is a purely linear variable, the constraint can't be written as a SOC */
2182 	         if( !SCIPisZero(scip, bp[i]) )
2183 	            goto CLEANUP;
2184 	
2185 	         bp[i] = 0.0;
2186 	         eigvals[i] = 0.0;
2187 	      }
2188 	      else if( eigvals[i] > 0.0 )
2189 	         npos++;
2190 	      else
2191 	         nneg++;
2192 	   }
2193 	
2194 	   /* a proper SOC constraint needs at least 2 variables */
2195 	   if( npos + nneg < 2 )
2196 	      goto CLEANUP;
2197 	
2198 	   /* determine whether rhs or lhs of cons is potentially SOC, if any */
2199 	   rhsissoc = (nneg == 1 && SCIPgetExprNLocksPosNonlinear(expr) > 0);
2200 	   lhsissoc = (npos == 1 && SCIPgetExprNLocksNegNonlinear(expr) > 0);
2201 	
2202 	   if( rhsissoc || lhsissoc )
2203 	   {
2204 	      if( conslhs == SCIP_INVALID || consrhs == SCIP_INVALID ) /*lint !e777*/
2205 	      {
2206 	         SCIP_CALL( SCIPevalExprActivity(scip, expr) );
2207 	         expractivity = SCIPexprGetActivity(expr);
2208 	         lhs = (conslhs == SCIP_INVALID ? expractivity.inf : conslhs); /*lint !e777*/
2209 	         rhs = (consrhs == SCIP_INVALID ? expractivity.sup : consrhs); /*lint !e777*/
2210 	      }
2211 	      else
2212 	      {
2213 	         lhs = conslhs;
2214 	         rhs = consrhs;
2215 	      }
2216 	   }
2217 	   else
2218 	   {
2219 	      /* if none of the sides is potentially SOC, stop */
2220 	      goto CLEANUP;
2221 	   }
2222 	
2223 	   /* @TODO: what do we do if both sides are possible? */
2224 	   if( !rhsissoc )
2225 	   {
2226 	      assert(lhsissoc);
2227 	
2228 	      /* lhs is potentially SOC, change signs */
2229 	      lhsconstant = lhs - constant;  /*lint !e644*/
2230 	
2231 	      for( i = 0; i < nvars; ++i )
2232 	      {
2233 	         eigvals[i] = -eigvals[i];
2234 	         bp[i] = -bp[i];
2235 	      }
2236 	      *enforcebelow = FALSE; /* enforce lhs <= expr */
2237 	   }
2238 	   else
2239 	   {
2240 	      lhsconstant = constant - rhs;  /*lint !e644*/
2241 	      *enforcebelow = TRUE; /* enforce expr <= rhs */
2242 	   }
2243 	
2244 	   /* initialize remaining datastructures for nonlinear handler */
2245 	   SCIP_CALL( SCIPallocBufferArray(scip, &offsets, npos + nneg + 1) );
2246 	   SCIP_CALL( SCIPallocBufferArray(scip, &transcoefs, ntranscoefs) );
2247 	   SCIP_CALL( SCIPallocBufferArray(scip, &transcoefsidx, ntranscoefs) );
2248 	   SCIP_CALL( SCIPallocBufferArray(scip, &termbegins, npos + nneg + 2) );
2249 	
2250 	   /* try to fill the nlhdlrexprdata (at this point, it can still fail) */
2251 	   SCIP_CALL( tryFillNlhdlrExprDataQuad(scip, occurringexprs, eigvecmatrix, eigvals, bp, nvars, termbegins, transcoefs,
2252 	         transcoefsidx, offsets, &lhsconstant, &nterms, success) );
2253 	
2254 	   if( !(*success) )
2255 	      goto CLEANUP;
2256 	
2257 	   assert(0 < nterms && nterms <= npos + nneg + 1);
2258 	   assert(ntranscoefs == termbegins[nterms]);
2259 	
2260 	   /*
2261 	    * at this point, the expression passed all checks and is SOC-representable
2262 	    */
2263 	
2264 	   /* register all requests for auxiliary variables */
2265 	   for( i = 0; i < nvars; ++i )
2266 	   {
2267 	      SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, occurringexprs[i], TRUE, FALSE, FALSE, FALSE) );
2268 	   }
2269 	
2270 	#ifdef SCIP_DEBUG
2271 	   SCIPdebugMsg(scip, "found SOC structure for expression %p\n%f <= ", (void*)expr, lhs);
2272 	   SCIPprintExpr(scip, expr, NULL);
2273 	   SCIPinfoMessage(scip, NULL, "<= %f\n", rhs);
2274 	#endif
2275 	
2276 	   /* finally, create and store nonlinear handler expression data */
2277 	   SCIP_CALL( createNlhdlrExprData(scip, occurringexprs, offsets, transcoefs, transcoefsidx, termbegins, nvars, nterms,
2278 	            nlhdlrexprdata) );
2279 	   assert(*nlhdlrexprdata != NULL);
2280 	
2281 	CLEANUP:
2282 	   SCIPfreeBufferArrayNull(scip, &termbegins);
2283 	   SCIPfreeBufferArrayNull(scip, &transcoefsidx);
2284 	   SCIPfreeBufferArrayNull(scip, &transcoefs);
2285 	   SCIPfreeBufferArrayNull(scip, &offsets);
2286 	   SCIPfreeBufferArrayNull(scip, &bp);
2287 	   SCIPfreeBufferArray(scip, &eigvals);
2288 	   SCIPfreeBufferArray(scip, &lincoefs);
2289 	   SCIPfreeBufferArray(scip, &eigvecmatrix);
2290 	   SCIPfreeBufferArray(scip, &occurringexprs);
2291 	   SCIPhashmapFree(&expr2idx);
2292 	
2293 	   return SCIP_OKAY;
2294 	}
2295 	
2296 	/** helper method to detect SOC structures
2297 	 *
2298 	 * The detection runs in 3 steps:
2299 	 *  1. check if expression is a norm of the form \f$\sqrt{\sum_i (\text{sqrcoef}_i\, \text{expr}_i^2 + \text{lincoef}_i\, \text{expr}_i) + c}\f$
2300 	 *  which can be transformed to the form \f$\sqrt{\sum_i (\text{coef}_i \text{expr}_i + \text{const}_i)^2 + c^*}\f$ with \f$c^* \geq 0\f$.\n
2301 	 *    -> this results in the SOC     expr &le; auxvar(expr)
2302 	 *
2303 	 *    TODO we should generalize and check for sqrt(positive-semidefinite-quadratic)
2304 	 *
2305 	 *  2. check if expression represents a quadratic function of one of the following forms (all coefs > 0)
2306 	 *     1. \f$(\sum_i   \text{coef}_i \text{expr}_i^2) - \text{coef}_k \text{expr}_k^2 \leq \text{RHS}\f$ or
2307 	 *     2. \f$(\sum_i - \text{coef}_i \text{expr}_i^2) + \text{coef}_k \text{expr}_k^2 \geq \text{LHS}\f$ or
2308 	 *     3. \f$(\sum_i   \text{coef}_i \text{expr}_i^2) - \text{coef}_k \text{expr}_k \text{expr}_l \leq \text{RHS}\f$ or
2309 	 *     4. \f$(\sum_i - \text{coef}_i \text{expr}_i^2) + \text{coef}_k \text{expr}_k \text{expr}_l \geq \text{LHS}\f$,
2310 	 *
2311 	 *     where RHS &ge; 0 or LHS &le; 0, respectively. For LHS and RHS we use the constraint sides if it is a root expr
2312 	 *     and the bounds of the auxiliary variable otherwise.
2313 	 *     The last two cases are called hyperbolic or rotated second order cone.\n
2314 	 *     -> this results in the SOC \f$\sqrt{(\sum_i \text{coef}_i \text{expr}_i^2) - \text{RHS}} \leq \sqrt{\text{coef}_k} \text{expr}_k\f$
2315 	 *                            or  \f$\sqrt{4(\sum_i \text{coef}_i \text{expr}_i^2) - 4\text{RHS} + (\text{expr}_k - \text{expr}_l)^2)} \leq \text{expr}_k + \text{expr}_l\f$.
2316 	 *                            (analogously for the LHS cases)
2317 	 *
2318 	 *  3. check if expression represents a quadratic inequality of the form \f$f(x) = x^TAx + b^Tx + c \leq 0\f$ such that \f$f(x)\f$
2319 	 *  has exactly one negative eigenvalue plus some extra conditions, see detectSocQuadraticComplex().
2320 	 *
2321 	 *  Note that step 3 is only performed if parameter `compeigenvalues` is set to TRUE.
2322 	 */
2323 	static
2324 	SCIP_RETCODE detectSOC(
2325 	   SCIP*                 scip,               /**< SCIP data structure */
2326 	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< nonlinear handler data */
2327 	   SCIP_EXPR*            expr,               /**< expression */
2328 	   SCIP_Real             conslhs,            /**< lhs of the constraint that the expression defines (or SCIP_INVALID) */
2329 	   SCIP_Real             consrhs,            /**< rhs of the constraint that the expression defines (or SCIP_INVALID) */
2330 	   SCIP_NLHDLREXPRDATA** nlhdlrexprdata,     /**< pointer to store nonlinear handler expression data */
2331 	   SCIP_Bool*            enforcebelow,       /**< pointer to store whether we enforce <= (TRUE) or >= (FALSE); only
2332 	                                              *   valid when success is TRUE */
2333 	   SCIP_Bool*            success             /**< pointer to store whether SOC structure has been detected */
2334 	   )
2335 	{
2336 	   assert(expr != NULL);
2337 	   assert(nlhdlrdata != NULL);
2338 	   assert(nlhdlrexprdata != NULL);
2339 	   assert(success != NULL);
2340 	
2341 	   *success = FALSE;
2342 	
2343 	   /* check whether expression is given as norm as described in case 1 above: if we have a constraint
2344 	    * sqrt(sum x_i^2) <= constant, then it might be better not to handle this here; thus, we only call detectSocNorm
2345 	    * when the expr is _not_ the root of a constraint
2346 	    */
2347 	   if( conslhs == SCIP_INVALID && consrhs == SCIP_INVALID ) /*lint !e777*/
2348 	   {
2349 	      SCIP_CALL( detectSocNorm(scip, expr, nlhdlrexprdata, success) );
2350 	      *enforcebelow = *success;
2351 	   }
2352 	
2353 	   if( !(*success) )
2354 	   {
2355 	      /* check whether expression is a simple soc-respresentable quadratic expression as described in case 2 above */
2356 	      SCIP_CALL( detectSocQuadraticSimple(scip, expr, conslhs, consrhs, nlhdlrexprdata, enforcebelow, success) );
2357 	   }
2358 	
2359 	   if( !(*success) && nlhdlrdata->compeigenvalues )
2360 	   {
2361 	      /* check whether expression is a more complex soc-respresentable quadratic expression as described in case 3 */
2362 	      SCIP_CALL( detectSocQuadraticComplex(scip, expr, conslhs, consrhs, nlhdlrexprdata, enforcebelow, success) );
2363 	   }
2364 	
2365 	   return SCIP_OKAY;
2366 	}
2367 	
2368 	/*
2369 	 * Callback methods of nonlinear handler
2370 	 */
2371 	
2372 	/** nonlinear handler copy callback */
2373 	static
2374 	SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrSoc)
2375 	{ /*lint --e{715}*/
2376 	   assert(targetscip != NULL);
2377 	   assert(sourcenlhdlr != NULL);
2378 	   assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), NLHDLR_NAME) == 0);
2379 	
2380 	   SCIP_CALL( SCIPincludeNlhdlrSoc(targetscip) );
2381 	
2382 	   return SCIP_OKAY;
2383 	}
2384 	
2385 	/** callback to free data of handler */
2386 	static
2387 	SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrFreehdlrdataSoc)
2388 	{ /*lint --e{715}*/
2389 	   assert(nlhdlrdata != NULL);
2390 	
2391 	   SCIPfreeBlockMemory(scip, nlhdlrdata);
2392 	
2393 	   return SCIP_OKAY;
2394 	}
2395 	
2396 	/** callback to free expression specific data */
2397 	static
2398 	SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataSoc)
2399 	{  /*lint --e{715}*/
2400 	   assert(*nlhdlrexprdata != NULL);
2401 	
2402 	   SCIP_CALL( freeNlhdlrExprData(scip, nlhdlrexprdata) );
2403 	
2404 	   return SCIP_OKAY;
2405 	}
2406 	
2407 	
2408 	/** callback to be called in initialization */
2409 	#if 0
2410 	static
2411 	SCIP_DECL_NLHDLRINIT(nlhdlrInitSoc)
2412 	{  /*lint --e{715}*/
2413 	   SCIPerrorMessage("method of soc nonlinear handler not implemented yet\n");
2414 	   SCIPABORT(); /*lint --e{527}*/
2415 	
2416 	   return SCIP_OKAY;
2417 	}
2418 	#else
2419 	#define nlhdlrInitSoc NULL
2420 	#endif
2421 	
2422 	
2423 	/** callback to be called in deinitialization */
2424 	#if 0
2425 	static
2426 	SCIP_DECL_NLHDLREXIT(nlhdlrExitSoc)
2427 	{  /*lint --e{715}*/
2428 	   SCIPerrorMessage("method of soc nonlinear handler not implemented yet\n");
2429 	   SCIPABORT(); /*lint --e{527}*/
2430 	
2431 	   return SCIP_OKAY;
2432 	}
2433 	#else
2434 	#define nlhdlrExitSoc NULL
2435 	#endif
2436 	
2437 	
2438 	/** callback to detect structure in expression tree */
2439 	static
2440 	SCIP_DECL_NLHDLRDETECT(nlhdlrDetectSoc)
2441 	{ /*lint --e{715}*/
2442 	   SCIP_Real conslhs;
2443 	   SCIP_Real consrhs;
2444 	   SCIP_Bool enforcebelow;
2445 	   SCIP_Bool success;
2446 	   SCIP_NLHDLRDATA* nlhdlrdata;
2447 	
2448 	   assert(expr != NULL);
2449 	
2450 	   /* don't try if no sepa is required
2451 	    * TODO implement some bound strengthening
2452 	    */
2453 	   if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABOTH) == SCIP_NLHDLR_METHOD_SEPABOTH )
2454 	      return SCIP_OKAY;
2455 	
2456 	   assert(SCIPgetExprNAuxvarUsesNonlinear(expr) > 0);  /* since some sepa is required, there should have been demand for it */
2457 	
2458 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
2459 	   assert(nlhdlrdata != NULL);
2460 	
2461 	   conslhs = (cons == NULL ? SCIP_INVALID : SCIPgetLhsNonlinear(cons));
2462 	   consrhs = (cons == NULL ? SCIP_INVALID : SCIPgetRhsNonlinear(cons));
2463 	
2464 	   SCIP_CALL( detectSOC(scip, nlhdlrdata, expr, conslhs, consrhs, nlhdlrexprdata, &enforcebelow, &success) );
2465 	
2466 	   if( !success )
2467 	      return SCIP_OKAY;
2468 	
2469 	   /* inform what we can do */
2470 	   *participating = enforcebelow ? SCIP_NLHDLR_METHOD_SEPABELOW : SCIP_NLHDLR_METHOD_SEPAABOVE;
2471 	
2472 	   /* if we have been successful on sqrt(...) <= auxvar, then we enforce
2473 	    * otherwise, expr is quadratic and we separate for expr <= ub(auxvar) only
2474 	    * in that case, we enforce only if expr is the root of a constraint, since then replacing auxvar by up(auxvar) does not relax anything (auxvar <= ub(auxvar) is the only constraint on auxvar)
2475 	    */
2476 	   if( (SCIPisExprPower(scip, expr) && SCIPgetExponentExprPow(expr) == 0.5) || (cons != NULL) )
2477 	      *enforcing |= *participating;
2478 	
2479 	   return SCIP_OKAY;
2480 	}
2481 	
2482 	
2483 	/** auxiliary evaluation callback of nonlinear handler
2484 	 * @todo: remember if we are in the original variables and avoid reevaluating
2485 	 */
2486 	static
2487 	SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxSoc)
2488 	{ /*lint --e{715}*/
2489 	   int i;
2490 	
2491 	   assert(nlhdlrexprdata != NULL);
2492 	   assert(nlhdlrexprdata->vars != NULL);
2493 	   assert(nlhdlrexprdata->transcoefs != NULL);
2494 	   assert(nlhdlrexprdata->transcoefsidx != NULL);
2495 	   assert(nlhdlrexprdata->nterms > 1);
2496 	
2497 	   /* if the original expression is a norm, evaluate w.r.t. the auxiliary variables */
2498 	   if( SCIPisExprPower(scip, expr) )
2499 	   {
2500 	      assert(SCIPgetExponentExprPow(expr) == 0.5);
2501 	
2502 	      updateVarVals(scip, nlhdlrexprdata, sol, FALSE);
2503 	
2504 	      /* compute sum_i coef_i expr_i^2 */
2505 	      *auxvalue = 0.0;
2506 	      for( i = 0; i < nlhdlrexprdata->nterms - 1; ++i )
2507 	      {
2508 	         SCIP_Real termval;
2509 	
2510 	         termval = evalSingleTerm(scip, nlhdlrexprdata, i);
2511 	         *auxvalue += SQR(termval);
2512 	      }
2513 	
2514 	      assert(*auxvalue >= 0.0);
2515 	
2516 	      /* compute SQRT(sum_i coef_i expr_i^2) */
2517 	      *auxvalue = SQRT(*auxvalue);
2518 	   }
2519 	   /* otherwise, evaluate the original quadratic expression w.r.t. the created auxvars of the children */
2520 	   else
2521 	   {
2522 	      SCIP_EXPR** children;
2523 	      SCIP_Real* childcoefs;
2524 	      int nchildren;
2525 	
2526 	      assert(SCIPisExprSum(scip, expr));
2527 	
2528 	      children = SCIPexprGetChildren(expr);
2529 	      childcoefs = SCIPgetCoefsExprSum(expr);
2530 	      nchildren = SCIPexprGetNChildren(expr);
2531 	
2532 	      *auxvalue = SCIPgetConstantExprSum(expr);
2533 	
2534 	      for( i = 0; i < nchildren; ++i )
2535 	      {
2536 	         if( SCIPisExprPower(scip, children[i]) )
2537 	         {
2538 	            SCIP_VAR* argauxvar;
2539 	            SCIP_Real solval;
2540 	
2541 	            assert(SCIPgetExponentExprPow(children[i]) == 2.0);
2542 	
2543 	            argauxvar = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(children[i])[0]);
2544 	            assert(argauxvar != NULL);
2545 	
2546 	            solval = SCIPgetSolVal(scip, sol, argauxvar);
2547 	            *auxvalue += childcoefs[i] * SQR( solval );
2548 	         }
2549 	         else if( SCIPisExprProduct(scip, children[i]) )
2550 	         {
2551 	            SCIP_VAR* argauxvar1;
2552 	            SCIP_VAR* argauxvar2;
2553 	
2554 	            assert(SCIPexprGetNChildren(children[i]) == 2);
2555 	
2556 	            argauxvar1 = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(children[i])[0]);
2557 	            argauxvar2 = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(children[i])[1]);
2558 	            assert(argauxvar1 != NULL);
2559 	            assert(argauxvar2 != NULL);
2560 	
2561 	            *auxvalue += childcoefs[i] * SCIPgetSolVal(scip, sol, argauxvar1) * SCIPgetSolVal(scip, sol, argauxvar2);
2562 	         }
2563 	         else
2564 	         {
2565 	            SCIP_VAR* argauxvar;
2566 	
2567 	            argauxvar = SCIPgetExprAuxVarNonlinear(children[i]);
2568 	            assert(argauxvar != NULL);
2569 	
2570 	            *auxvalue += childcoefs[i] * SCIPgetSolVal(scip, sol, argauxvar);
2571 	         }
2572 	      }
2573 	   }
2574 	
2575 	   return SCIP_OKAY;
2576 	}
2577 	
2578 	
2579 	/** separation deinitialization method of a nonlinear handler (called during CONSINITLP) */
2580 	static
2581 	SCIP_DECL_NLHDLRINITSEPA(nlhdlrInitSepaSoc)
2582 	{ /*lint --e{715}*/
2583 	   SCIP_ROWPREP* rowprep;
2584 	
2585 	   assert(conshdlr != NULL);
2586 	   assert(expr != NULL);
2587 	   assert(nlhdlrexprdata != NULL);
2588 	
2589 	   /* already needed for debug solution */
2590 	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &nlhdlrexprdata->varvals, nlhdlrexprdata->nvars) );
2591 	
2592 	   /* if we have 3 or more terms in lhs create variable and row for disaggregation */
2593 	   if( nlhdlrexprdata->nterms > 3 )
2594 	   {
2595 	      /* create variables for cone disaggregation */
2596 	      SCIP_CALL( createDisaggrVars(scip, expr, nlhdlrexprdata) );
2597 	
2598 	#ifdef WITH_DEBUG_SOLUTION
2599 	      if( SCIPdebugIsMainscip(scip) )
2600 	      {
2601 	         SCIP_Real lhsval;
2602 	         SCIP_Real rhsval;
2603 	         SCIP_Real disvarval;
2604 	         int ndisvars;
2605 	         int nterms;
2606 	         int i;
2607 	
2608 	         for( i = 0; i < nlhdlrexprdata->nvars; ++i )
2609 	         {
2610 	            SCIP_CALL( SCIPdebugGetSolVal(scip, SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[i]), &nlhdlrexprdata->varvals[i]) );
2611 	         }
2612 	
2613 	         /*  the debug solution value of the disaggregation variables is set to
2614 	          *      (v_i^T x + beta_i)^2 / (v_{n+1}^T x + beta_{n+1})
2615 	          *  if (v_{n+1}^T x + beta_{n+1}) is different from 0.
2616 	          *  Otherwise, the debug solution value is set to 0.
2617 	          */
2618 	
2619 	         nterms = nlhdlrexprdata->nterms;
2620 	
2621 	         /* find value of rhs */
2622 	         rhsval = evalSingleTerm(scip, nlhdlrexprdata, nterms - 1);
2623 	
2624 	         /* set value of disaggregation vars */
2625 	         ndisvars = nlhdlrexprdata->nterms - 1;
2626 	
2627 	         if( SCIPisZero(scip, rhsval) )
2628 	         {
2629 	            for( i = 0; i < ndisvars; ++i )
2630 	            {
2631 	               SCIP_CALL( SCIPdebugAddSolVal(scip, nlhdlrexprdata->disvars[i], 0.0) );
2632 	            }
2633 	         }
2634 	         else
2635 	         {
2636 	            /* set value for each disaggregation variable corresponding to quadratic term */
2637 	            for( i = 0; i < ndisvars; ++i )
2638 	            {
2639 	               lhsval = evalSingleTerm(scip, nlhdlrexprdata, i);
2640 	
2641 	               disvarval = SQR(lhsval) / rhsval;
2642 	
2643 	               SCIP_CALL( SCIPdebugAddSolVal(scip, nlhdlrexprdata->disvars[i], disvarval) );
2644 	            }
2645 	         }
2646 	      }
2647 	#endif
2648 	
2649 	      /* create the disaggregation row and store it in nlhdlrexprdata */
2650 	      SCIP_CALL( createDisaggrRow(scip, conshdlr, expr, nlhdlrexprdata) );
2651 	   }
2652 	
2653 	#ifdef SCIP_DEBUG
2654 	   SCIPdebugMsg(scip, "initlp for \n");
2655 	   printNlhdlrExprData(scip, nlhdlrexprdata);
2656 	#endif
2657 	
2658 	   /* add some initial cuts on well-selected coordinates */
2659 	   if( nlhdlrexprdata->nterms == 2 )
2660 	   {
2661 	      /* we have |v_1^T x + \beta_1| \leq v_2^T x + \beta_2
2662 	       *
2663 	       * we should linearize at points where the first term is -1 or 1, so we can take
2664 	       *
2665 	       * x = v_1 / ||v_1||^2 (+/-1 - beta_1)
2666 	       */
2667 	      SCIP_Real plusminus1;
2668 	      SCIP_Real norm;
2669 	      int i;
2670 	
2671 	      /* calculate ||v_1||^2 */
2672 	      norm = 0.0;
2673 	      for( i = nlhdlrexprdata->termbegins[0]; i < nlhdlrexprdata->termbegins[1]; ++i )
2674 	         norm += SQR(nlhdlrexprdata->transcoefs[i]);
2675 	      assert(norm > 0.0);
2676 	
2677 	      BMSclearMemoryArray(nlhdlrexprdata->varvals, nlhdlrexprdata->nvars);
2678 	
2679 	      for( plusminus1 = -1.0; plusminus1 <= 1.0; plusminus1 += 2.0 )
2680 	      {
2681 	         /* set x = v_1 / ||v_1||^2 (plusminus1 - beta_1) */
2682 	         for( i = nlhdlrexprdata->termbegins[0]; i < nlhdlrexprdata->termbegins[1]; ++i )
2683 	            nlhdlrexprdata->varvals[nlhdlrexprdata->transcoefsidx[i]] = nlhdlrexprdata->transcoefs[i] / norm * (plusminus1 - nlhdlrexprdata->offsets[0]);
2684 	         assert(SCIPisEQ(scip, evalSingleTerm(scip, nlhdlrexprdata, 0), plusminus1));
2685 	
2686 	         /* compute gradient cut */
2687 	         SCIP_CALL( generateCutSolSOC(scip, &rowprep, expr, cons, nlhdlrexprdata, -SCIPinfinity(scip), evalSingleTerm(scip, nlhdlrexprdata, 1)) );
2688 	
2689 	         if( rowprep != NULL )
2690 	         {
2691 	            SCIP_Bool success = FALSE;
2692 	
2693 	            SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, NULL, SCIPgetHugeValue(scip), &success) );
2694 	            if( success )
2695 	            {
2696 	               SCIP_ROW* cut;
2697 	               SCIP_CALL( SCIPgetRowprepRowCons(scip, &cut, rowprep, cons) );
2698 	               SCIP_CALL( SCIPaddRow(scip, cut, FALSE, infeasible) );
2699 	               SCIP_CALL( SCIPreleaseRow(scip, &cut) );
2700 	            }
2701 	
2702 	            SCIPfreeRowprep(scip, &rowprep);
2703 	         }
2704 	
2705 	         if( *infeasible )
2706 	            return SCIP_OKAY;
2707 	      }
2708 	   }
2709 	   else if( nlhdlrexprdata->nterms == 3 && nlhdlrexprdata->termbegins[0] != nlhdlrexprdata->termbegins[1] )
2710 	   {
2711 	      /* we have \sqrt{ (v_1^T x + \beta_1)^2 + (v_2^T x + \beta_2)^2 } \leq v_3^T x + \beta_3
2712 	       * with v_1 != 0
2713 	       *
2714 	       * we should linearize at points where the first and second term are (-1,0), (1,1), (1,-1), i.e.,
2715 	       *
2716 	       * v_1^T x + \beta_1 = -1   1   1
2717 	       * v_2^T x + \beta_2 =  0   1  -1
2718 	       *
2719 	       * Let i be the index of the first nonzero element in v_1.
2720 	       * Let j != i be an index of v_2, to be determined.
2721 	       * Assume all other entries of x will be set to 0.
2722 	       * Then we have
2723 	       *   (v_1)_i x_i + (v_1)_j x_j = c   with  c = -1 - beta_1
2724 	       *   (v_2)_i x_i + (v_2)_j x_j = d   with  d =  0 - beta_2
2725 	       *
2726 	       * Since (v_1)_i != 0, this gives
2727 	       *   x_i = 1/(v_1)_i (c - (v_1)_j x_j)
2728 	       * Substituting in the 2nd equation, we get
2729 	       *    (v_2)_i/(v_1)_i (c - (v_1)_j x_j) + (v_2)_j x_j = d
2730 	       * -> ((v_2)_j - (v_2)_i (v_1)_j / (v_1)_i) x_j = d - (v_2)_i/(v_1)_i c
2731 	       * Now find j such that (v_2)_j - (v_2)_i (v_1)_j / (v_1)_i != 0.
2732 	       *
2733 	       * If v_2 = 0, then linearize only for first term being -1 or 1 and don't care about value of second term.
2734 	       * We then set j arbitrary, x_i = 1/(v_1)_i c, other coordinates of x = zero.
2735 	       */
2736 	      static const SCIP_Real refpoints[3][2] = { {-1.0, 0.0}, {1.0, 1.0}, {1.0, -1.0} };
2737 	      SCIP_Real v1i, v1j = 0.0;
2738 	      SCIP_Real v2i, v2j = 0.0;
2739 	      SCIP_Bool v2zero;
2740 	      int i;
2741 	      int j = -1;
2742 	      int k;
2743 	      int pos;
2744 	
2745 	      i = nlhdlrexprdata->transcoefsidx[nlhdlrexprdata->termbegins[0]];
2746 	      v1i = nlhdlrexprdata->transcoefs[nlhdlrexprdata->termbegins[0]];
2747 	      assert(v1i != 0.0);
2748 	
2749 	      v2zero = nlhdlrexprdata->termbegins[1] == nlhdlrexprdata->termbegins[2];
2750 	
2751 	      /* get (v_2)_i; as it is a sparse vector, we need to search for i in transcoefsidx */
2752 	      v2i = 0.0;
2753 	      if( !v2zero && SCIPsortedvecFindInt(nlhdlrexprdata->transcoefsidx + nlhdlrexprdata->termbegins[1], i, nlhdlrexprdata->termbegins[2] - nlhdlrexprdata->termbegins[1], &pos) )
2754 	      {
2755 	         assert(nlhdlrexprdata->transcoefsidx[nlhdlrexprdata->termbegins[1] + pos] == i);
2756 	         v2i = nlhdlrexprdata->transcoefs[nlhdlrexprdata->termbegins[1] + pos];
2757 	      }
2758 	
2759 	      /* find a j, for now look only into indices with (v_2)_j != 0 */
2760 	      for( k = nlhdlrexprdata->termbegins[1]; k < nlhdlrexprdata->termbegins[2]; ++k )
2761 	      {
2762 	         /* check whether transcoefsidx[k] could be a good j */
2763 	
2764 	         if( nlhdlrexprdata->transcoefsidx[k] == i ) /* i == j is not good */
2765 	            continue;
2766 	
2767 	         /* get (v_1)_j; as it is a sparse vector, we need to search for j in transcoefsidx */
2768 	         v1j = 0.0;
2769 	         if( SCIPsortedvecFindInt(nlhdlrexprdata->transcoefsidx + nlhdlrexprdata->termbegins[0], nlhdlrexprdata->transcoefsidx[k], nlhdlrexprdata->termbegins[1] - nlhdlrexprdata->termbegins[0], &pos) )
2770 	         {
2771 	            assert(nlhdlrexprdata->transcoefsidx[nlhdlrexprdata->termbegins[0] + pos] == nlhdlrexprdata->transcoefsidx[k]);
2772 	            v1j = nlhdlrexprdata->transcoefs[nlhdlrexprdata->termbegins[0] + pos];
2773 	         }
2774 	
2775 	         v2j = nlhdlrexprdata->transcoefs[k];
2776 	
2777 	         if( SCIPisZero(scip, v2j - v2i * v1j / v1i) )  /* (v_2)_j - (v_2)_i (v_1)_j / (v_1)_i = 0 is also not good */
2778 	            continue;
2779 	
2780 	         j = nlhdlrexprdata->transcoefsidx[k];
2781 	         break;
2782 	      }
2783 	
2784 	      if( v2zero )
2785 	      {
2786 	         j = 0;
2787 	         v1j = 0.0;
2788 	         v2j = 0.0;
2789 	      }
2790 	
2791 	      if( j != -1 )
2792 	      {
2793 	         SCIP_Real c, d;
2794 	         int point;
2795 	
2796 	         BMSclearMemoryArray(nlhdlrexprdata->varvals, nlhdlrexprdata->nvars);
2797 	
2798 	         for( point = 0; point < (v2zero ? 2 : 3); ++point )
2799 	         {
2800 	            c = refpoints[point][0] - nlhdlrexprdata->offsets[0];
2801 	
2802 	            if( !v2zero )
2803 	            {
2804 	               /* set x_j and x_i */
2805 	               d = refpoints[point][1] - nlhdlrexprdata->offsets[1];
2806 	               nlhdlrexprdata->varvals[j] = (d - v2i/v1i*c) / (v2j - v2i * v1j / v1i);
2807 	               nlhdlrexprdata->varvals[i] = (c - v1j * nlhdlrexprdata->varvals[j]) / v1i;
2808 	
2809 	               SCIPdebugMsg(scip, "<%s>(%d) = %g, <%s>(%d) = %g\n",
2810 	                  SCIPvarGetName(SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[i])), i, nlhdlrexprdata->varvals[i],
2811 	                  SCIPvarGetName(SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[j])), j, nlhdlrexprdata->varvals[j]);
2812 	            }
2813 	            else
2814 	            {
2815 	               /* set x_i */
2816 	               nlhdlrexprdata->varvals[i] = c / v1i;
2817 	
2818 	               SCIPdebugMsg(scip, "<%s>(%d) = %g\n",
2819 	                  SCIPvarGetName(SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[i])), i, nlhdlrexprdata->varvals[i]);
2820 	            }
2821 	
2822 	            assert(SCIPisEQ(scip, evalSingleTerm(scip, nlhdlrexprdata, 0), refpoints[point][0]));
2823 	            assert(v2zero || SCIPisEQ(scip, evalSingleTerm(scip, nlhdlrexprdata, 1), refpoints[point][1]));
2824 	
2825 	            /* compute gradient cut */
2826 	            SCIP_CALL( generateCutSolSOC(scip, &rowprep, expr, cons, nlhdlrexprdata, -SCIPinfinity(scip), evalSingleTerm(scip, nlhdlrexprdata, 2)) );
2827 	
2828 	            if( rowprep != NULL )
2829 	            {
2830 	               SCIP_Bool success = FALSE;
2831 	
2832 	               SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, NULL, SCIPgetHugeValue(scip), &success) );
2833 	               if( success )
2834 	               {
2835 	                  SCIP_ROW* cut;
2836 	                  SCIP_CALL( SCIPgetRowprepRowCons(scip, &cut, rowprep, cons) );
2837 	                  SCIP_CALL( SCIPaddRow(scip, cut, FALSE, infeasible) );
2838 	                  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
2839 	               }
2840 	
2841 	               SCIPfreeRowprep(scip, &rowprep);
2842 	            }
2843 	
2844 	            if( *infeasible )
2845 	               return SCIP_OKAY;
2846 	         }
2847 	      }
2848 	   }
2849 	   else if( nlhdlrexprdata->nterms == 3 )
2850 	   {
2851 	      /* we have \sqrt{ \beta_1^2 + (v_2^T x + \beta_2)^2 } \leq v_3^T x + \beta_3
2852 	       * with v_2 != 0
2853 	       *
2854 	       * we should linearize at points where the second term is -1 or 1
2855 	       *
2856 	       * set x = v_2 / ||v_2||^2 (+/-1 - beta_2)
2857 	       */
2858 	      SCIP_Real plusminus1;
2859 	      SCIP_Real norm;
2860 	      int i;
2861 	
2862 	      /* calculate ||v_2||^2 */
2863 	      norm = 0.0;
2864 	      for( i = nlhdlrexprdata->termbegins[1]; i < nlhdlrexprdata->termbegins[2]; ++i )
2865 	         norm += SQR(nlhdlrexprdata->transcoefs[i]);
2866 	      assert(norm > 0.0);
2867 	
2868 	      BMSclearMemoryArray(nlhdlrexprdata->varvals, nlhdlrexprdata->nvars);
2869 	
2870 	      for( plusminus1 = -1.0; plusminus1 <= 1.0; plusminus1 += 2.0 )
2871 	      {
2872 	         /* set x = v_2 / ||v_2||^2 (plusminus1 - beta_2) */
2873 	         for( i = nlhdlrexprdata->termbegins[1]; i < nlhdlrexprdata->termbegins[2]; ++i )
2874 	            nlhdlrexprdata->varvals[nlhdlrexprdata->transcoefsidx[i]] = nlhdlrexprdata->transcoefs[i] / norm * (plusminus1 - nlhdlrexprdata->offsets[1]);
2875 	         assert(SCIPisEQ(scip, evalSingleTerm(scip, nlhdlrexprdata, 1), plusminus1));
2876 	
2877 	         /* compute gradient cut */
2878 	         SCIP_CALL( generateCutSolSOC(scip, &rowprep, expr, cons, nlhdlrexprdata, -SCIPinfinity(scip), evalSingleTerm(scip, nlhdlrexprdata, 2)) );
2879 	
2880 	         if( rowprep != NULL )
2881 	         {
2882 	            SCIP_Bool success = FALSE;
2883 	
2884 	            SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, NULL, SCIPgetHugeValue(scip), &success) );
2885 	            if( success )
2886 	            {
2887 	               SCIP_ROW* cut;
2888 	               SCIP_CALL( SCIPgetRowprepRowCons(scip, &cut, rowprep, cons) );
2889 	               SCIP_CALL( SCIPaddRow(scip, cut, FALSE, infeasible) );
2890 	               SCIP_CALL( SCIPreleaseRow(scip, &cut) );
2891 	            }
2892 	
2893 	            SCIPfreeRowprep(scip, &rowprep);
2894 	         }
2895 	
2896 	         if( *infeasible )
2897 	            return SCIP_OKAY;
2898 	      }
2899 	   }
2900 	   else
2901 	   {
2902 	      /* generate gradient cuts for the small rotated cones
2903 	       *  \f[
2904 	       *    \sqrt{4(v_k^T x + \beta_k)^2 + (v_n^T x + \beta_n - y_k)^2} - v_n^T x - \beta_n - y_k \leq 0.
2905 	       *  \f]
2906 	       *
2907 	       * we should linearize again at points where the first and second term (inside sqr) are (-1/2,0), (1/2,1), (1/2,-1).
2908 	       * Since we have y_k, we can achieve this more easily here via
2909 	       *   x = v_k/||v_k||^2 (+/-0.5 - beta_k)
2910 	       *   y_k = v_n^T x + beta_n + 0/1/-1
2911 	       *
2912 	       * If v_k = 0, then we use x = 0 and linearize for second term being 1 and -1 only
2913 	       */
2914 	      static const SCIP_Real refpoints[3][2] = { {-0.5, 0.0}, {0.5, 1.0}, {0.5, -1.0} };
2915 	      SCIP_Real rhsval;
2916 	      SCIP_Real norm;
2917 	      int point;
2918 	      int i;
2919 	      int k;
2920 	      SCIP_Bool vkzero;
2921 	
2922 	      /* add disaggregation row to LP */
2923 	      SCIP_CALL( SCIPaddRow(scip, nlhdlrexprdata->disrow, FALSE, infeasible) );
2924 	
2925 	      if( *infeasible )
2926 	         return SCIP_OKAY;
2927 	
2928 	      for( k = 0; k < nlhdlrexprdata->nterms - 1; ++k )
2929 	      {
2930 	         vkzero = nlhdlrexprdata->termbegins[k+1] == nlhdlrexprdata->termbegins[k];
2931 	         assert(!vkzero || nlhdlrexprdata->offsets[k] != 0.0);
2932 	
2933 	         /* calculate ||v_k||^2 */
2934 	         norm = 0.0;
2935 	         for( i = nlhdlrexprdata->termbegins[k]; i < nlhdlrexprdata->termbegins[k+1]; ++i )
2936 	            norm += SQR(nlhdlrexprdata->transcoefs[i]);
2937 	         assert(vkzero || norm > 0.0);
2938 	
2939 	         BMSclearMemoryArray(nlhdlrexprdata->varvals, nlhdlrexprdata->nvars);
2940 	
2941 	         for( point = vkzero ? 1 : 0; point < 3; ++point )
2942 	         {
2943 	            /* set x = v_k / ||v_k||^2 (refpoints[point][0] - beta_k) / 2 */
2944 	            for( i = nlhdlrexprdata->termbegins[k]; i < nlhdlrexprdata->termbegins[k+1]; ++i )
2945 	               nlhdlrexprdata->varvals[nlhdlrexprdata->transcoefsidx[i]] = nlhdlrexprdata->transcoefs[i] / norm * (refpoints[point][0] - nlhdlrexprdata->offsets[k]);  /*lint !e795*/
2946 	            assert(vkzero || SCIPisEQ(scip, evalSingleTerm(scip, nlhdlrexprdata, k), refpoints[point][0]));
2947 	
2948 	            /* set y_k = v_n^T x + beta_n + 0/1/-1 */
2949 	            rhsval = evalSingleTerm(scip, nlhdlrexprdata, nlhdlrexprdata->nterms - 1);
2950 	            nlhdlrexprdata->disvarvals[k] = rhsval + refpoints[point][1];
2951 	
2952 	            /* compute gradient cut */
2953 	            SCIP_CALL( generateCutSolDisagg(scip, &rowprep, expr, cons, nlhdlrexprdata, k, -SCIPinfinity(scip), rhsval) );
2954 	
2955 	            if( rowprep != NULL )
2956 	            {
2957 	               SCIP_Bool success = FALSE;
2958 	
2959 	               SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, NULL, SCIPgetHugeValue(scip), &success) );
2960 	               if( success )
2961 	               {
2962 	                  SCIP_ROW* cut;
2963 	                  SCIP_CALL( SCIPgetRowprepRowCons(scip, &cut, rowprep, cons) );
2964 	                  SCIP_CALL( SCIPaddRow(scip, cut, FALSE, infeasible) );
2965 	                  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
2966 	               }
2967 	
2968 	               SCIPfreeRowprep(scip, &rowprep);
2969 	            }
2970 	
2971 	            if( *infeasible )
2972 	               return SCIP_OKAY;
2973 	         }
2974 	      }
2975 	   }
2976 	
2977 	   return SCIP_OKAY;
2978 	}
2979 	
2980 	
2981 	/** separation deinitialization method of a nonlinear handler (called during CONSEXITSOL) */
2982 	static
2983 	SCIP_DECL_NLHDLREXITSEPA(nlhdlrExitSepaSoc)
2984 	{ /*lint --e{715}*/
2985 	   assert(nlhdlrexprdata != NULL);
2986 	
2987 	   /* free disaggreagation row */
2988 	   if( nlhdlrexprdata->disrow != NULL )
2989 	   {
2990 	      SCIP_CALL( SCIPreleaseRow(scip, &nlhdlrexprdata->disrow) );
2991 	   }
2992 	
2993 	   SCIPfreeBlockMemoryArray(scip, &nlhdlrexprdata->varvals, nlhdlrexprdata->nvars);
2994 	
2995 	   return SCIP_OKAY;
2996 	}
2997 	
2998 	
2999 	/** nonlinear handler separation callback */
3000 	static
3001 	SCIP_DECL_NLHDLRENFO(nlhdlrEnfoSoc)
3002 	{ /*lint --e{715}*/
3003 	   SCIP_NLHDLRDATA* nlhdlrdata;
3004 	   SCIP_Real rhsval;
3005 	   int ndisaggrs;
3006 	   int k;
3007 	   SCIP_Bool infeasible;
3008 	
3009 	   assert(nlhdlrexprdata != NULL);
3010 	   assert(nlhdlrexprdata->nterms < 4 || nlhdlrexprdata->disrow != NULL);
3011 	   assert(nlhdlrexprdata->nterms > 1);
3012 	
3013 	   *result = SCIP_DIDNOTFIND;
3014 	
3015 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
3016 	   assert(nlhdlrdata != NULL);
3017 	
3018 	   /* update varvals
3019 	    * set variables close to integer to integer, in particular when close to zero
3020 	    * for simple soc's (no large v_i, no offsets), variables close to zero would give coefficients close to zero in the cut,
3021 	    * which the cut cleanup may have problems to relax (and we end up with local or much relaxed cuts)
3022 	    * also when close to other integers, rounding now may prevent some relaxation in cut cleanup
3023 	    */
3024 	   updateVarVals(scip, nlhdlrexprdata, sol, TRUE);
3025 	
3026 	   rhsval = evalSingleTerm(scip, nlhdlrexprdata, nlhdlrexprdata->nterms - 1);
3027 	
3028 	   /* if there are three or two terms just compute gradient cut */
3029 	   if( nlhdlrexprdata->nterms < 4 )
3030 	   {
3031 	      SCIP_ROWPREP* rowprep;
3032 	
3033 	      /* compute gradient cut */
3034 	      SCIP_CALL( generateCutSolSOC(scip, &rowprep, expr, cons, nlhdlrexprdata, SCIPgetLPFeastol(scip), rhsval) );
3035 	
3036 	      if( rowprep != NULL )
3037 	      {
3038 	         SCIP_CALL( addCut(scip, nlhdlrdata, rowprep, sol, cons, allowweakcuts, result) );
3039 	
3040 	         SCIPfreeRowprep(scip, &rowprep);
3041 	      }
3042 	      else
3043 	      {
3044 	         SCIPdebugMsg(scip, "failed to generate cut for SOC\n");
3045 	      }
3046 	
3047 	      return SCIP_OKAY;
3048 	   }
3049 	
3050 	   ndisaggrs = nlhdlrexprdata->nterms - 1;
3051 	
3052 	   /* check whether the aggregation row is in the LP */
3053 	   if( !SCIProwIsInLP(nlhdlrexprdata->disrow) && -SCIPgetRowSolFeasibility(scip, nlhdlrexprdata->disrow, sol) > SCIPgetLPFeastol(scip) )
3054 	   {
3055 	      SCIP_CALL( SCIPaddRow(scip, nlhdlrexprdata->disrow, TRUE, &infeasible) );
3056 	      SCIPdebugMsg(scip, "added disaggregation row to LP, cutoff=%u\n", infeasible);
3057 	
3058 	      if( infeasible )
3059 	      {
3060 	         *result = SCIP_CUTOFF;
3061 	         return SCIP_OKAY;
3062 	      }
3063 	
3064 	      *result = SCIP_SEPARATED;
3065 	   }
3066 	
3067 	   for( k = 0; k < ndisaggrs && *result != SCIP_CUTOFF; ++k )
3068 	   {
3069 	      SCIP_ROWPREP* rowprep;
3070 	
3071 	      /* compute gradient cut */
3072 	      SCIP_CALL( generateCutSolDisagg(scip, &rowprep, expr, cons, nlhdlrexprdata, k, SCIPgetLPFeastol(scip), rhsval) );
3073 	
3074 	      if( rowprep != NULL )
3075 	      {
3076 	         SCIP_CALL( addCut(scip, nlhdlrdata, rowprep, sol, cons, allowweakcuts, result) );
3077 	
3078 	         SCIPfreeRowprep(scip, &rowprep);
3079 	      }
3080 	   }
3081 	
3082 	   return SCIP_OKAY;
3083 	}
3084 	
3085 	static
3086 	SCIP_DECL_NLHDLRSOLLINEARIZE(nlhdlrSollinearizeSoc)
3087 	{ /*lint --e{715}*/
3088 	   SCIP_NLHDLRDATA* nlhdlrdata;
3089 	   SCIP_Real rhsval;
3090 	   int k;
3091 	
3092 	   assert(sol != NULL);
3093 	   assert(nlhdlrexprdata != NULL);
3094 	   assert(nlhdlrexprdata->nterms < 4 || nlhdlrexprdata->disrow != NULL);
3095 	   assert(nlhdlrexprdata->nterms > 1);
3096 	
3097 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
3098 	   assert(nlhdlrdata != NULL);
3099 	
3100 	   /* update varvals */
3101 	   updateVarVals(scip, nlhdlrexprdata, sol, TRUE);
3102 	
3103 	   rhsval = evalSingleTerm(scip, nlhdlrexprdata, nlhdlrexprdata->nterms - 1);
3104 	
3105 	   /* if there are three or two terms just compute gradient cut */
3106 	   if( nlhdlrexprdata->nterms < 4 )
3107 	   {
3108 	      SCIP_ROWPREP* rowprep;
3109 	
3110 	      /* compute gradient cut */
3111 	      SCIP_CALL( generateCutSolSOC(scip, &rowprep, expr, cons, nlhdlrexprdata, -SCIPinfinity(scip), rhsval) );
3112 	
3113 	      if( rowprep != NULL )
3114 	      {
3115 	         SCIP_CALL( addCutPool(scip, nlhdlrdata, rowprep, sol, cons) );
3116 	
3117 	         SCIPfreeRowprep(scip, &rowprep);
3118 	      }
3119 	
3120 	      return SCIP_OKAY;
3121 	   }
3122 	
3123 	   for( k = 0; k < nlhdlrexprdata->nterms - 1; ++k )
3124 	   {
3125 	      SCIP_ROWPREP* rowprep;
3126 	
3127 	      /* compute gradient cut */
3128 	      SCIP_CALL( generateCutSolDisagg(scip, &rowprep, expr, cons, nlhdlrexprdata, k, -SCIPinfinity(scip), rhsval) );
3129 	
3130 	      if( rowprep != NULL )
3131 	      {
3132 	         SCIP_CALL( addCutPool(scip, nlhdlrdata, rowprep, sol, cons) );
3133 	
3134 	         SCIPfreeRowprep(scip, &rowprep);
3135 	      }
3136 	   }
3137 	
3138 	   return SCIP_OKAY;
3139 	}
3140 	
3141 	/*
3142 	 * nonlinear handler specific interface methods
3143 	 */
3144 	
3145 	/** includes SOC nonlinear handler in nonlinear constraint handler */
3146 	SCIP_RETCODE SCIPincludeNlhdlrSoc(
3147 	   SCIP*                 scip                /**< SCIP data structure */
3148 	   )
3149 	{
3150 	   SCIP_NLHDLRDATA* nlhdlrdata;
3151 	   SCIP_NLHDLR* nlhdlr;
3152 	
3153 	   assert(scip != NULL);
3154 	
3155 	   /* create nonlinear handler data */
3156 	   SCIP_CALL( SCIPallocClearBlockMemory(scip, &nlhdlrdata) );
3157 	
3158 	   SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY, NLHDLR_ENFOPRIORITY, nlhdlrDetectSoc, nlhdlrEvalauxSoc, nlhdlrdata) );
3159 	   assert(nlhdlr != NULL);
3160 	
3161 	   SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrSoc);
3162 	   SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrFreehdlrdataSoc);
3163 	   SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeExprDataSoc);
3164 	   SCIPnlhdlrSetInitExit(nlhdlr, nlhdlrInitSoc, nlhdlrExitSoc);
3165 	   SCIPnlhdlrSetSepa(nlhdlr, nlhdlrInitSepaSoc, nlhdlrEnfoSoc, NULL, nlhdlrExitSepaSoc);
3166 	   SCIPnlhdlrSetSollinearize(nlhdlr, nlhdlrSollinearizeSoc);
3167 	
3168 	   /* add soc nlhdlr parameters */
3169 	   /* TODO should we get rid of this and use separating/mineffiacy(root) instead, which is 1e-4? */
3170 	   SCIP_CALL( SCIPaddRealParam(scip, "nlhdlr/" NLHDLR_NAME "/mincutefficacy",
3171 	         "Minimum efficacy which a cut needs in order to be added.",
3172 	         &nlhdlrdata->mincutefficacy, FALSE, DEFAULT_MINCUTEFFICACY, 0.0, SCIPinfinity(scip), NULL, NULL) );
3173 	
3174 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/compeigenvalues",
3175 	         "Should Eigenvalue computations be done to detect complex cases in quadratic constraints?",
3176 	         &nlhdlrdata->compeigenvalues, FALSE, DEFAULT_COMPEIGENVALUES, NULL, NULL) );
3177 	
3178 	   return SCIP_OKAY;
3179 	}
3180 	
3181 	/** checks whether constraint is SOC representable in original variables and returns the SOC representation
3182 	 *
3183 	 * The SOC representation has the form:
3184 	 * \f$\sqrt{\sum_{i=1}^{n} (v_i^T x + \beta_i)^2} - v_{n+1}^T x - \beta_{n+1} \lessgtr 0\f$,
3185 	 * where \f$n+1 = \text{nterms}\f$ and the inequality type is given by sidetype (`SCIP_SIDETYPE_RIGHT` if inequality
3186 	 * is \f$\leq\f$, `SCIP_SIDETYPE_LEFT` if \f$\geq\f$).
3187 	 *
3188 	 * For each term (i.e. for each \f$i\f$ in the above notation as well as \f$n+1\f$), the constant \f$\beta_i\f$ is given by the
3189 	 * corresponding element `offsets[i-1]` and `termbegins[i-1]` is the starting position of the term in arrays
3190 	 * `transcoefs` and `transcoefsidx`. The overall number of nonzeros is `termbegins[nterms]`.
3191 	 *
3192 	 * Arrays `transcoefs` and `transcoefsidx` have size `termbegins[nterms]` and define the linear expressions \f$v_i^T x\f$
3193 	 * for each term. For a term \f$i\f$ in the above notation, the nonzeroes are given by elements
3194 	 * `termbegins[i-1]...termbegins[i]` of `transcoefs` and `transcoefsidx`. There may be no nonzeroes for some term (i.e.,
3195 	 * constant terms are possible). `transcoefs` contains the coefficients \f$v_i\f$ and `transcoefsidx` contains positions of
3196 	 * variables in the `vars` array.
3197 	 *
3198 	 * The `vars` array has size `nvars` and contains \f$x\f$ variables; each variable is included at most once.
3199 	 *
3200 	 * The arrays should be freed by calling SCIPfreeSOCArraysNonlinear().
3201 	 *
3202 	 * This function uses the methods that are used in the detection algorithm of the SOC nonlinear handler.
3203 	 */
3204 	SCIP_RETCODE SCIPisSOCNonlinear(
3205 	   SCIP*                 scip,               /**< SCIP data structure */
3206 	   SCIP_CONS*            cons,               /**< nonlinear constraint */
3207 	   SCIP_Bool             compeigenvalues,    /**< whether eigenvalues should be computed to detect complex cases */
3208 	   SCIP_Bool*            success,            /**< pointer to store whether SOC structure has been detected */
3209 	   SCIP_SIDETYPE*        sidetype,           /**< pointer to store which side of cons is SOC representable; only
3210 	                                              *   valid when success is TRUE */
3211 	   SCIP_VAR***           vars,               /**< variables (x) that appear on both sides; no duplicates are allowed */
3212 	   SCIP_Real**           offsets,            /**< offsets of both sides (beta_i) */
3213 	   SCIP_Real**           transcoefs,         /**< non-zeros of linear transformation vectors (v_i) */
3214 	   int**                 transcoefsidx,      /**< mapping of transformation coefficients to variable indices in vars */
3215 	   int**                 termbegins,         /**< starting indices of transcoefs for each term */
3216 	   int*                  nvars,              /**< total number of variables appearing (i.e. size of vars) */
3217 	   int*                  nterms              /**< number of summands in the SQRT +1 for RHS (n+1) */
3218 	   )
3219 	{
3220 	   SCIP_NLHDLRDATA nlhdlrdata;
3221 	   SCIP_NLHDLREXPRDATA *nlhdlrexprdata;
3222 	   SCIP_Real conslhs;
3223 	   SCIP_Real consrhs;
3224 	   SCIP_EXPR* expr;
3225 	   SCIP_Bool enforcebelow;
3226 	   int i;
3227 	
3228 	   assert(cons != NULL);
3229 	
3230 	   expr = SCIPgetExprNonlinear(cons);
3231 	   assert(expr != NULL);
3232 	
3233 	   nlhdlrdata.mincutefficacy = 0.0;
3234 	   nlhdlrdata.compeigenvalues = compeigenvalues;
3235 	
3236 	   conslhs = SCIPgetLhsNonlinear(cons);
3237 	   consrhs = SCIPgetRhsNonlinear(cons);
3238 	
3239 	   SCIP_CALL( detectSOC(scip, &nlhdlrdata, expr, conslhs, consrhs, &nlhdlrexprdata, &enforcebelow, success) );
3240 	
3241 	   /* the constraint must be SOC representable in original variables */
3242 	   if( *success )
3243 	   {
3244 	      assert(nlhdlrexprdata != NULL);
3245 	
3246 	      for( i = 0; i < nlhdlrexprdata->nvars; ++i )
3247 	      {
3248 	         if( !SCIPisExprVar(scip, nlhdlrexprdata->vars[i]) )
3249 	         {
3250 	            *success = FALSE;
3251 	            break;
3252 	         }
3253 	      }
3254 	   }
3255 	
3256 	   if( *success )
3257 	   {
3258 	      *sidetype = enforcebelow ? SCIP_SIDETYPE_RIGHT : SCIP_SIDETYPE_LEFT;
3259 	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, vars, nlhdlrexprdata->nvars) );
3260 	
3261 	      for( i = 0; i < nlhdlrexprdata->nvars; ++i )
3262 	      {
3263 	         (*vars)[i] = SCIPgetVarExprVar(nlhdlrexprdata->vars[i]);
3264 	         assert((*vars)[i] != NULL);
3265 	      }
3266 	      SCIPfreeBlockMemoryArray(scip, &nlhdlrexprdata->vars, nlhdlrexprdata->nvars);
3267 	      *offsets = nlhdlrexprdata->offsets;
3268 	      *transcoefs = nlhdlrexprdata->transcoefs;
3269 	      *transcoefsidx = nlhdlrexprdata->transcoefsidx;
3270 	      *termbegins = nlhdlrexprdata->termbegins;
3271 	      *nvars = nlhdlrexprdata->nvars;
3272 	      *nterms = nlhdlrexprdata->nterms;
3273 	      SCIPfreeBlockMemory(scip, &nlhdlrexprdata);
3274 	   }
3275 	   else
3276 	   {
3277 	      if( nlhdlrexprdata != NULL )
3278 	      {
3279 	         SCIP_CALL( freeNlhdlrExprData(scip, &nlhdlrexprdata) );
3280 	      }
3281 	      *vars = NULL;
3282 	      *offsets = NULL;
3283 	      *transcoefs = NULL;
3284 	      *transcoefsidx = NULL;
3285 	      *termbegins = NULL;
3286 	      *nvars = 0;
3287 	      *nterms = 0;
3288 	   }
3289 	
3290 	   return SCIP_OKAY;
3291 	}
3292 	
3293 	/** frees arrays created by SCIPisSOCNonlinear() */
3294 	void SCIPfreeSOCArraysNonlinear(
3295 	   SCIP*                 scip,               /**< SCIP data structure */
3296 	   SCIP_VAR***           vars,               /**< variables that appear on both sides (x) */
3297 	   SCIP_Real**           offsets,            /**< offsets of both sides (beta_i) */
3298 	   SCIP_Real**           transcoefs,         /**< non-zeros of linear transformation vectors (v_i) */
3299 	   int**                 transcoefsidx,      /**< mapping of transformation coefficients to variable indices in vars */
3300 	   int**                 termbegins,         /**< starting indices of transcoefs for each term */
3301 	   int                   nvars,              /**< total number of variables appearing */
3302 	   int                   nterms              /**< number of summands in the SQRT +1 for RHS (n+1) */
3303 	   )
3304 	{
3305 	   int ntranscoefs;
3306 	
3307 	   if( nvars == 0 )
3308 	      return;
3309 	
3310 	   assert(vars != NULL);
3311 	   assert(offsets != NULL);
3312 	   assert(transcoefs != NULL);
3313 	   assert(transcoefsidx != NULL);
3314 	   assert(termbegins != NULL);
3315 	
3316 	   ntranscoefs = (*termbegins)[nterms];
3317 	
3318 	   SCIPfreeBlockMemoryArray(scip, termbegins, nterms + 1);
3319 	   SCIPfreeBlockMemoryArray(scip, transcoefsidx, ntranscoefs);
3320 	   SCIPfreeBlockMemoryArray(scip, transcoefs, ntranscoefs);
3321 	   SCIPfreeBlockMemoryArray(scip, offsets, nterms);
3322 	   SCIPfreeBlockMemoryArray(scip, vars, nvars);
3323 	}
3324