Scippy

SCIP

Solving Constraint Integer Programs

nlhdlr_soc.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (c) 2002-2024 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 */
600  (void) SCIPsnprintf(SCIProwprepGetName(*rowprep), SCIP_MAXSTRLEN, "soc%d_%p_%" SCIP_LONGINT_FORMAT, 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_%" SCIP_LONGINT_FORMAT, (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  if( branchcandonly )
3016  return SCIP_OKAY;
3017 
3018  nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
3019  assert(nlhdlrdata != NULL);
3020 
3021  /* update varvals
3022  * set variables close to integer to integer, in particular when close to zero
3023  * for simple soc's (no large v_i, no offsets), variables close to zero would give coefficients close to zero in the cut,
3024  * which the cut cleanup may have problems to relax (and we end up with local or much relaxed cuts)
3025  * also when close to other integers, rounding now may prevent some relaxation in cut cleanup
3026  */
3027  updateVarVals(scip, nlhdlrexprdata, sol, TRUE);
3028 
3029  rhsval = evalSingleTerm(scip, nlhdlrexprdata, nlhdlrexprdata->nterms - 1);
3030 
3031  /* if there are three or two terms just compute gradient cut */
3032  if( nlhdlrexprdata->nterms < 4 )
3033  {
3034  SCIP_ROWPREP* rowprep;
3035 
3036  /* compute gradient cut */
3037  SCIP_CALL( generateCutSolSOC(scip, &rowprep, expr, cons, nlhdlrexprdata, SCIPgetLPFeastol(scip), rhsval) );
3038 
3039  if( rowprep != NULL )
3040  {
3041  SCIP_CALL( addCut(scip, nlhdlrdata, rowprep, sol, cons, allowweakcuts, result) );
3042 
3043  SCIPfreeRowprep(scip, &rowprep);
3044  }
3045  else
3046  {
3047  SCIPdebugMsg(scip, "failed to generate cut for SOC\n");
3048  }
3049 
3050  return SCIP_OKAY;
3051  }
3052 
3053  ndisaggrs = nlhdlrexprdata->nterms - 1;
3054 
3055  /* check whether the aggregation row is in the LP */
3056  if( !SCIProwIsInLP(nlhdlrexprdata->disrow) && -SCIPgetRowSolFeasibility(scip, nlhdlrexprdata->disrow, sol) > SCIPgetLPFeastol(scip) )
3057  {
3058  SCIP_CALL( SCIPaddRow(scip, nlhdlrexprdata->disrow, TRUE, &infeasible) );
3059  SCIPdebugMsg(scip, "added disaggregation row to LP, cutoff=%u\n", infeasible);
3060 
3061  if( infeasible )
3062  {
3063  *result = SCIP_CUTOFF;
3064  return SCIP_OKAY;
3065  }
3066 
3067  *result = SCIP_SEPARATED;
3068  }
3069 
3070  for( k = 0; k < ndisaggrs && *result != SCIP_CUTOFF; ++k )
3071  {
3072  SCIP_ROWPREP* rowprep;
3073 
3074  /* compute gradient cut */
3075  SCIP_CALL( generateCutSolDisagg(scip, &rowprep, expr, cons, nlhdlrexprdata, k, SCIPgetLPFeastol(scip), rhsval) );
3076 
3077  if( rowprep != NULL )
3078  {
3079  SCIP_CALL( addCut(scip, nlhdlrdata, rowprep, sol, cons, allowweakcuts, result) );
3080 
3081  SCIPfreeRowprep(scip, &rowprep);
3082  }
3083  }
3084 
3085  return SCIP_OKAY;
3086 }
3087 
3088 static
3089 SCIP_DECL_NLHDLRSOLLINEARIZE(nlhdlrSollinearizeSoc)
3090 { /*lint --e{715}*/
3091  SCIP_NLHDLRDATA* nlhdlrdata;
3092  SCIP_Real rhsval;
3093  int k;
3094 
3095  assert(sol != NULL);
3096  assert(nlhdlrexprdata != NULL);
3097  assert(nlhdlrexprdata->nterms < 4 || nlhdlrexprdata->disrow != NULL);
3098  assert(nlhdlrexprdata->nterms > 1);
3099 
3100  nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
3101  assert(nlhdlrdata != NULL);
3102 
3103  /* update varvals */
3104  updateVarVals(scip, nlhdlrexprdata, sol, TRUE);
3105 
3106  rhsval = evalSingleTerm(scip, nlhdlrexprdata, nlhdlrexprdata->nterms - 1);
3107 
3108  /* if there are three or two terms just compute gradient cut */
3109  if( nlhdlrexprdata->nterms < 4 )
3110  {
3111  SCIP_ROWPREP* rowprep;
3112 
3113  /* compute gradient cut */
3114  SCIP_CALL( generateCutSolSOC(scip, &rowprep, expr, cons, nlhdlrexprdata, -SCIPinfinity(scip), rhsval) );
3115 
3116  if( rowprep != NULL )
3117  {
3118  SCIP_CALL( addCutPool(scip, nlhdlrdata, rowprep, sol, cons) );
3119 
3120  SCIPfreeRowprep(scip, &rowprep);
3121  }
3122 
3123  return SCIP_OKAY;
3124  }
3125 
3126  for( k = 0; k < nlhdlrexprdata->nterms - 1; ++k )
3127  {
3128  SCIP_ROWPREP* rowprep;
3129 
3130  /* compute gradient cut */
3131  SCIP_CALL( generateCutSolDisagg(scip, &rowprep, expr, cons, nlhdlrexprdata, k, -SCIPinfinity(scip), rhsval) );
3132 
3133  if( rowprep != NULL )
3134  {
3135  SCIP_CALL( addCutPool(scip, nlhdlrdata, rowprep, sol, cons) );
3136 
3137  SCIPfreeRowprep(scip, &rowprep);
3138  }
3139  }
3140 
3141  return SCIP_OKAY;
3142 }
3143 
3144 /*
3145  * nonlinear handler specific interface methods
3146  */
3147 
3148 /** includes SOC nonlinear handler in nonlinear constraint handler */
3149 SCIP_RETCODE SCIPincludeNlhdlrSoc(
3150  SCIP* scip /**< SCIP data structure */
3151  )
3152 {
3153  SCIP_NLHDLRDATA* nlhdlrdata;
3154  SCIP_NLHDLR* nlhdlr;
3155 
3156  assert(scip != NULL);
3157 
3158  /* create nonlinear handler data */
3159  SCIP_CALL( SCIPallocClearBlockMemory(scip, &nlhdlrdata) );
3160 
3161  SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY, NLHDLR_ENFOPRIORITY, nlhdlrDetectSoc, nlhdlrEvalauxSoc, nlhdlrdata) );
3162  assert(nlhdlr != NULL);
3163 
3164  SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrSoc);
3165  SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrFreehdlrdataSoc);
3166  SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeExprDataSoc);
3167  SCIPnlhdlrSetInitExit(nlhdlr, nlhdlrInitSoc, nlhdlrExitSoc);
3168  SCIPnlhdlrSetSepa(nlhdlr, nlhdlrInitSepaSoc, nlhdlrEnfoSoc, NULL, nlhdlrExitSepaSoc);
3169  SCIPnlhdlrSetSollinearize(nlhdlr, nlhdlrSollinearizeSoc);
3170 
3171  /* add soc nlhdlr parameters */
3172  /* TODO should we get rid of this and use separating/mineffiacy(root) instead, which is 1e-4? */
3173  SCIP_CALL( SCIPaddRealParam(scip, "nlhdlr/" NLHDLR_NAME "/mincutefficacy",
3174  "Minimum efficacy which a cut needs in order to be added.",
3175  &nlhdlrdata->mincutefficacy, FALSE, DEFAULT_MINCUTEFFICACY, 0.0, SCIPinfinity(scip), NULL, NULL) );
3176 
3177  SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/compeigenvalues",
3178  "Should Eigenvalue computations be done to detect complex cases in quadratic constraints?",
3179  &nlhdlrdata->compeigenvalues, FALSE, DEFAULT_COMPEIGENVALUES, NULL, NULL) );
3180 
3181  return SCIP_OKAY;
3182 }
3183 
3184 /** checks whether constraint is SOC representable in original variables and returns the SOC representation
3185  *
3186  * The SOC representation has the form:
3187  * \f$\sqrt{\sum_{i=1}^{n} (v_i^T x + \beta_i)^2} - v_{n+1}^T x - \beta_{n+1} \lessgtr 0\f$,
3188  * where \f$n+1 = \text{nterms}\f$ and the inequality type is given by sidetype (`SCIP_SIDETYPE_RIGHT` if inequality
3189  * is \f$\leq\f$, `SCIP_SIDETYPE_LEFT` if \f$\geq\f$).
3190  *
3191  * 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
3192  * corresponding element `offsets[i-1]` and `termbegins[i-1]` is the starting position of the term in arrays
3193  * `transcoefs` and `transcoefsidx`. The overall number of nonzeros is `termbegins[nterms]`.
3194  *
3195  * Arrays `transcoefs` and `transcoefsidx` have size `termbegins[nterms]` and define the linear expressions \f$v_i^T x\f$
3196  * for each term. For a term \f$i\f$ in the above notation, the nonzeroes are given by elements
3197  * `termbegins[i-1]...termbegins[i]` of `transcoefs` and `transcoefsidx`. There may be no nonzeroes for some term (i.e.,
3198  * constant terms are possible). `transcoefs` contains the coefficients \f$v_i\f$ and `transcoefsidx` contains positions of
3199  * variables in the `vars` array.
3200  *
3201  * The `vars` array has size `nvars` and contains \f$x\f$ variables; each variable is included at most once.
3202  *
3203  * The arrays should be freed by calling SCIPfreeSOCArraysNonlinear().
3204  *
3205  * This function uses the methods that are used in the detection algorithm of the SOC nonlinear handler.
3206  */
3207 SCIP_RETCODE SCIPisSOCNonlinear(
3208  SCIP* scip, /**< SCIP data structure */
3209  SCIP_CONS* cons, /**< nonlinear constraint */
3210  SCIP_Bool compeigenvalues, /**< whether eigenvalues should be computed to detect complex cases */
3211  SCIP_Bool* success, /**< pointer to store whether SOC structure has been detected */
3212  SCIP_SIDETYPE* sidetype, /**< pointer to store which side of cons is SOC representable; only
3213  * valid when success is TRUE */
3214  SCIP_VAR*** vars, /**< variables (x) that appear on both sides; no duplicates are allowed */
3215  SCIP_Real** offsets, /**< offsets of both sides (beta_i) */
3216  SCIP_Real** transcoefs, /**< non-zeros of linear transformation vectors (v_i) */
3217  int** transcoefsidx, /**< mapping of transformation coefficients to variable indices in vars */
3218  int** termbegins, /**< starting indices of transcoefs for each term */
3219  int* nvars, /**< total number of variables appearing (i.e. size of vars) */
3220  int* nterms /**< number of summands in the SQRT +1 for RHS (n+1) */
3221  )
3222 {
3223  SCIP_NLHDLRDATA nlhdlrdata;
3224  SCIP_NLHDLREXPRDATA *nlhdlrexprdata;
3225  SCIP_Real conslhs;
3226  SCIP_Real consrhs;
3227  SCIP_EXPR* expr;
3228  SCIP_Bool enforcebelow;
3229  int i;
3230 
3231  assert(cons != NULL);
3232 
3233  expr = SCIPgetExprNonlinear(cons);
3234  assert(expr != NULL);
3235 
3236  nlhdlrdata.mincutefficacy = 0.0;
3237  nlhdlrdata.compeigenvalues = compeigenvalues;
3238 
3239  conslhs = SCIPgetLhsNonlinear(cons);
3240  consrhs = SCIPgetRhsNonlinear(cons);
3241 
3242  SCIP_CALL( detectSOC(scip, &nlhdlrdata, expr, conslhs, consrhs, &nlhdlrexprdata, &enforcebelow, success) );
3243 
3244  /* the constraint must be SOC representable in original variables */
3245  if( *success )
3246  {
3247  assert(nlhdlrexprdata != NULL);
3248 
3249  for( i = 0; i < nlhdlrexprdata->nvars; ++i )
3250  {
3251  if( !SCIPisExprVar(scip, nlhdlrexprdata->vars[i]) )
3252  {
3253  *success = FALSE;
3254  break;
3255  }
3256  }
3257  }
3258 
3259  if( *success )
3260  {
3261  *sidetype = enforcebelow ? SCIP_SIDETYPE_RIGHT : SCIP_SIDETYPE_LEFT;
3262  SCIP_CALL( SCIPallocBlockMemoryArray(scip, vars, nlhdlrexprdata->nvars) );
3263 
3264  for( i = 0; i < nlhdlrexprdata->nvars; ++i )
3265  {
3266  (*vars)[i] = SCIPgetVarExprVar(nlhdlrexprdata->vars[i]);
3267  assert((*vars)[i] != NULL);
3268  }
3269  SCIPfreeBlockMemoryArray(scip, &nlhdlrexprdata->vars, nlhdlrexprdata->nvars);
3270  *offsets = nlhdlrexprdata->offsets;
3271  *transcoefs = nlhdlrexprdata->transcoefs;
3272  *transcoefsidx = nlhdlrexprdata->transcoefsidx;
3273  *termbegins = nlhdlrexprdata->termbegins;
3274  *nvars = nlhdlrexprdata->nvars;
3275  *nterms = nlhdlrexprdata->nterms;
3276  SCIPfreeBlockMemory(scip, &nlhdlrexprdata);
3277  }
3278  else
3279  {
3280  if( nlhdlrexprdata != NULL )
3281  {
3282  SCIP_CALL( freeNlhdlrExprData(scip, &nlhdlrexprdata) );
3283  }
3284  *vars = NULL;
3285  *offsets = NULL;
3286  *transcoefs = NULL;
3287  *transcoefsidx = NULL;
3288  *termbegins = NULL;
3289  *nvars = 0;
3290  *nterms = 0;
3291  }
3292 
3293  return SCIP_OKAY;
3294 }
3295 
3296 /** frees arrays created by SCIPisSOCNonlinear() */
3297 void SCIPfreeSOCArraysNonlinear(
3298  SCIP* scip, /**< SCIP data structure */
3299  SCIP_VAR*** vars, /**< variables that appear on both sides (x) */
3300  SCIP_Real** offsets, /**< offsets of both sides (beta_i) */
3301  SCIP_Real** transcoefs, /**< non-zeros of linear transformation vectors (v_i) */
3302  int** transcoefsidx, /**< mapping of transformation coefficients to variable indices in vars */
3303  int** termbegins, /**< starting indices of transcoefs for each term */
3304  int nvars, /**< total number of variables appearing */
3305  int nterms /**< number of summands in the SQRT +1 for RHS (n+1) */
3306  )
3307 {
3308  int ntranscoefs;
3309 
3310  if( nvars == 0 )
3311  return;
3312 
3313  assert(vars != NULL);
3314  assert(offsets != NULL);
3315  assert(transcoefs != NULL);
3316  assert(transcoefsidx != NULL);
3317  assert(termbegins != NULL);
3318 
3319  ntranscoefs = (*termbegins)[nterms];
3320 
3321  SCIPfreeBlockMemoryArray(scip, termbegins, nterms + 1);
3322  SCIPfreeBlockMemoryArray(scip, transcoefsidx, ntranscoefs);
3323  SCIPfreeBlockMemoryArray(scip, transcoefs, ntranscoefs);
3324  SCIPfreeBlockMemoryArray(scip, offsets, nterms);
3325  SCIPfreeBlockMemoryArray(scip, vars, nvars);
3326 }
SCIP_RESULT
enum SCIP_Result SCIP_RESULT
Definition: type_result.h:61
SCIPfreeBlockMemoryArray
#define SCIPfreeBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:110
updateVarVals
static void updateVarVals(SCIP *scip, SCIP_NLHDLREXPRDATA *nlhdlrexprdata, SCIP_SOL *sol, SCIP_Bool roundtinyfrac)
Definition: nlhdlr_soc.c:412
nlhdlrInitSoc
#define nlhdlrInitSoc
Definition: nlhdlr_soc.c:2419
SCIP_DECL_NLHDLRFREEHDLRDATA
static SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrFreehdlrdataSoc)
Definition: nlhdlr_soc.c:2387
nterms
static volatile int nterms
Definition: interrupt.c:47
NULL
#define NULL
Definition: def.h:267
SCIPallocBlockMemoryArray
#define SCIPallocBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:93
SCIP_Interval
Definition: intervalarith.h:53
createDisaggrVars
static SCIP_RETCODE createDisaggrVars(SCIP *scip, SCIP_EXPR *expr, SCIP_NLHDLREXPRDATA *nlhdlrexprdata)
Definition: nlhdlr_soc.c:216
SCIPhashsetRemove
SCIP_RETCODE SCIPhashsetRemove(SCIP_HASHSET *hashset, void *element)
Definition: misc.c:3858
SCIPprintExpr
SCIP_RETCODE SCIPprintExpr(SCIP *scip, SCIP_EXPR *expr, FILE *file)
Definition: scip_expr.c:1486
SCIPlapackComputeEigenvalues
SCIP_RETCODE SCIPlapackComputeEigenvalues(BMS_BUFMEM *bufmem, SCIP_Bool geteigenvectors, int N, SCIP_Real *a, SCIP_Real *w)
Definition: lapack_calls.c:352
SCIPgetExprNAuxvarUsesNonlinear
unsigned int SCIPgetExprNAuxvarUsesNonlinear(SCIP_EXPR *expr)
Definition: cons_nonlinear.c:14378
SCIPevalExprActivity
SCIP_RETCODE SCIPevalExprActivity(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1717
SCIP_NLHDLR_METHOD_SEPABELOW
#define SCIP_NLHDLR_METHOD_SEPABELOW
Definition: type_nlhdlr.h:51
SCIPhashsetGetNElements
int SCIPhashsetGetNElements(SCIP_HASHSET *hashset)
Definition: misc.c:3992
SCIPexprGetNChildren
int SCIPexprGetNChildren(SCIP_EXPR *expr)
Definition: expr.c:3854
SCIPallocClearBufferArray
#define SCIPallocClearBufferArray(scip, ptr, num)
Definition: scip_mem.h:126
tryFillNlhdlrExprDataQuad
static SCIP_RETCODE tryFillNlhdlrExprDataQuad(SCIP *scip, SCIP_EXPR **occurringexprs, SCIP_Real *eigvecmatrix, SCIP_Real *eigvals, SCIP_Real *bp, int nvars, int *termbegins, SCIP_Real *transcoefs, int *transcoefsidx, SCIP_Real *offsets, SCIP_Real *lhsconstant, int *nterms, SCIP_Bool *success)
Definition: nlhdlr_soc.c:1071
SCIP_MAXSTRLEN
#define SCIP_MAXSTRLEN
Definition: def.h:288
SCIPnlhdlrSetFreeExprData
void SCIPnlhdlrSetFreeExprData(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRFREEEXPRDATA((*freeexprdata)))
Definition: nlhdlr.c:98
addCutPool
static SCIP_RETCODE addCutPool(SCIP *scip, SCIP_NLHDLRDATA *nlhdlrdata, SCIP_ROWPREP *rowprep, SCIP_SOL *sol, SCIP_CONS *cons)
Definition: nlhdlr_soc.c:831
SCIPcleanupRowprep2
SCIP_RETCODE SCIPcleanupRowprep2(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_SOL *sol, SCIP_Real maxcoefbound, SCIP_Bool *success)
Definition: misc_rowprep.c:1376
SCIPaddVarToRow
SCIP_RETCODE SCIPaddVarToRow(SCIP *scip, SCIP_ROW *row, SCIP_VAR *var, SCIP_Real val)
Definition: scip_lp.c:1701
SQR
#define SQR(x)
Definition: def.h:214
SCIPisPositive
SCIP_Bool SCIPisPositive(SCIP *scip, SCIP_Real val)
Definition: scip_numerics.c:570
SCIPgetConstantExprSum
SCIP_Real SCIPgetConstantExprSum(SCIP_EXPR *expr)
Definition: expr_sum.c:1567
addCut
static SCIP_RETCODE addCut(SCIP *scip, SCIP_NLHDLRDATA *nlhdlrdata, SCIP_ROWPREP *rowprep, SCIP_SOL *sol, SCIP_CONS *cons, SCIP_Bool allowweakcuts, SCIP_RESULT *result)
Definition: nlhdlr_soc.c:763
NLHDLR_DETECTPRIORITY
#define NLHDLR_DETECTPRIORITY
Definition: nlhdlr_soc.c:58
SCIPgetRhsNonlinear
SCIP_Real SCIPgetRhsNonlinear(SCIP_CONS *cons)
Definition: cons_nonlinear.c:13737
SCIP_HashSet
Definition: struct_misc.h:149
SCIP_Var
Definition: struct_var.h:207
SCIPreleaseVar
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip_var.c:1247
SCIPvarIsBinary
SCIP_Bool SCIPvarIsBinary(SCIP_VAR *var)
Definition: var.c:17600
SCIPregisterExprUsageNonlinear
SCIP_RETCODE SCIPregisterExprUsageNonlinear(SCIP *scip, SCIP_EXPR *expr, SCIP_Bool useauxvar, SCIP_Bool useactivityforprop, SCIP_Bool useactivityforsepabelow, SCIP_Bool useactivityforsepaabove)
Definition: cons_nonlinear.c:14401
FALSE
#define FALSE
Definition: def.h:94
SCIPhashmapCreate
SCIP_RETCODE SCIPhashmapCreate(SCIP_HASHMAP **hashmap, BMS_BLKMEM *blkmem, int mapsize)
Definition: misc.c:3074
SCIPnlhdlrSetFreeHdlrData
void SCIPnlhdlrSetFreeHdlrData(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRFREEHDLRDATA((*freehdlrdata)))
Definition: nlhdlr.c:87
SCIProwprepGetName
char * SCIProwprepGetName(SCIP_ROWPREP *rowprep)
Definition: misc_rowprep.c:689
SCIPinfinity
SCIP_Real SCIPinfinity(SCIP *scip)
Definition: scip_numerics.c:510
SCIPsnprintf
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition: misc.c:10877
SCIPisNegative
SCIP_Bool SCIPisNegative(SCIP *scip, SCIP_Real val)
Definition: scip_numerics.c:582
SCIPgetExponentExprPow
SCIP_Real SCIPgetExponentExprPow(SCIP_EXPR *expr)
Definition: expr_pow.c:3457
TRUE
#define TRUE
Definition: def.h:93
SCIPhashsetCreate
SCIP_RETCODE SCIPhashsetCreate(SCIP_HASHSET **hashset, BMS_BLKMEM *blkmem, int size)
Definition: misc.c:3759
SCIP_RETCODE
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:63
DEFAULT_COMPEIGENVALUES
#define DEFAULT_COMPEIGENVALUES
Definition: nlhdlr_soc.c:61
SCIPhashmapInsertInt
SCIP_RETCODE SCIPhashmapInsertInt(SCIP_HASHMAP *hashmap, void *origin, int image)
Definition: misc.c:3192
SCIPgetLPFeastol
SCIP_Real SCIPgetLPFeastol(SCIP *scip)
Definition: scip_lp.c:428
SCIPhashsetExists
SCIP_Bool SCIPhashsetExists(SCIP_HASHSET *hashset, void *element)
Definition: misc.c:3817
SCIPcreateVarBasic
SCIP_RETCODE SCIPcreateVarBasic(SCIP *scip, SCIP_VAR **var, const char *name, SCIP_Real lb, SCIP_Real ub, SCIP_Real obj, SCIP_VARTYPE vartype)
Definition: scip_var.c:194
SCIPexprGetActivity
SCIP_INTERVAL SCIPexprGetActivity(SCIP_EXPR *expr)
Definition: expr.c:4010
SCIPswapReals
void SCIPswapReals(SCIP_Real *value1, SCIP_Real *value2)
Definition: misc.c:10383
SCIPfreeBlockMemory
#define SCIPfreeBlockMemory(scip, ptr)
Definition: scip_mem.h:108
SCIPfreeSOCArraysNonlinear
void SCIPfreeSOCArraysNonlinear(SCIP *scip, SCIP_VAR ***vars, SCIP_Real **offsets, SCIP_Real **transcoefs, int **transcoefsidx, int **termbegins, int nvars, int nterms)
Definition: nlhdlr_soc.c:3297
SCIP_SIDETYPE_LEFT
Definition: type_lp.h:64
SCIP_SEPARATED
Definition: type_result.h:49
SCIPduplicateBufferArray
#define SCIPduplicateBufferArray(scip, ptr, source, num)
Definition: scip_mem.h:132
createNlhdlrExprData
static SCIP_RETCODE createNlhdlrExprData(SCIP *scip, SCIP_EXPR **vars, SCIP_Real *offsets, SCIP_Real *transcoefs, int *transcoefsidx, int *termbegins, int nvars, int nterms, SCIP_NLHDLREXPRDATA **nlhdlrexprdata)
Definition: nlhdlr_soc.c:335
SCIPisEQ
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:445
detectSocQuadraticComplex
static SCIP_RETCODE detectSocQuadraticComplex(SCIP *scip, SCIP_EXPR *expr, SCIP_Real conslhs, SCIP_Real consrhs, SCIP_NLHDLREXPRDATA **nlhdlrexprdata, SCIP_Bool *enforcebelow, SCIP_Bool *success)
Definition: nlhdlr_soc.c:2050
SCIP_VERBLEVEL_FULL
Definition: type_message.h:57
SCIPfreeBufferArray
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:136
SCIPallocBlockMemory
#define SCIPallocBlockMemory(scip, ptr)
Definition: scip_mem.h:89
expr_var.h
variable expression handler
SCIPensureRowprepSize
SCIP_RETCODE SCIPensureRowprepSize(SCIP *scip, SCIP_ROWPREP *rowprep, int size)
Definition: misc_rowprep.c:887
SCIPdebugMsg
#define SCIPdebugMsg
Definition: scip_message.h:78
detectSocQuadraticSimple
static SCIP_RETCODE detectSocQuadraticSimple(SCIP *scip, SCIP_EXPR *expr, SCIP_Real conslhs, SCIP_Real consrhs, SCIP_NLHDLREXPRDATA **nlhdlrexprdata, SCIP_Bool *enforcebelow, SCIP_Bool *success)
Definition: nlhdlr_soc.c:1582
SCIP_Nlhdlr
Definition: struct_nlhdlr.h:43
SCIPinfoMessage
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:208
SCIPisExprProduct
SCIP_Bool SCIPisExprProduct(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1464
SCIPgetRowprepRowCons
SCIP_RETCODE SCIPgetRowprepRowCons(SCIP *scip, SCIP_ROW **row, SCIP_ROWPREP *rowprep, SCIP_CONS *cons)
Definition: misc_rowprep.c:1678
SCIP_Sol
Definition: struct_sol.h:73
SCIP_DECL_NLHDLRSOLLINEARIZE
static SCIP_DECL_NLHDLRSOLLINEARIZE(nlhdlrSollinearizeSoc)
Definition: nlhdlr_soc.c:3089
SCIProwIsInLP
SCIP_Bool SCIProwIsInLP(SCIP_ROW *row)
Definition: lp.c:17523
SCIPhashmapExists
SCIP_Bool SCIPhashmapExists(SCIP_HASHMAP *hashmap, void *origin)
Definition: misc.c:3423
SCIPaddVarLocksType
SCIP_RETCODE SCIPaddVarLocksType(SCIP *scip, SCIP_VAR *var, SCIP_LOCKTYPE locktype, int nlocksdown, int nlocksup)
Definition: scip_var.c:4258
detectSocNorm
static SCIP_RETCODE detectSocNorm(SCIP *scip, SCIP_EXPR *expr, SCIP_NLHDLREXPRDATA **nlhdlrexprdata, SCIP_Bool *success)
Definition: nlhdlr_soc.c:1298
generateCutSolDisagg
static SCIP_RETCODE generateCutSolDisagg(SCIP *scip, SCIP_ROWPREP **rowprep, SCIP_EXPR *expr, SCIP_CONS *cons, SCIP_NLHDLREXPRDATA *nlhdlrexprdata, int disaggidx, SCIP_Real mincutviolation, SCIP_Real rhsval)
Definition: nlhdlr_soc.c:627
SCIPexprGetChildren
SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr)
Definition: expr.c:3864
DEFAULT_MINCUTEFFICACY
#define DEFAULT_MINCUTEFFICACY
Definition: nlhdlr_soc.c:60
SCIPduplicateBlockMemoryArray
#define SCIPduplicateBlockMemoryArray(scip, ptr, source, num)
Definition: scip_mem.h:105
SCIP_HashMap
Definition: struct_misc.h:137
SCIP_Interval::inf
SCIP_Real inf
Definition: intervalarith.h:55
SCIPbuffer
BMS_BUFMEM * SCIPbuffer(SCIP *scip)
Definition: scip_mem.c:72
NLHDLR_ENFOPRIORITY
#define NLHDLR_ENFOPRIORITY
Definition: nlhdlr_soc.c:59
SCIPmarkRowNotRemovableLocal
void SCIPmarkRowNotRemovableLocal(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1868
SCIPgetCoefsExprSum
SCIP_Real * SCIPgetCoefsExprSum(SCIP_EXPR *expr)
Definition: expr_sum.c:1552
SCIP_DECL_NLHDLRINITSEPA
static SCIP_DECL_NLHDLRINITSEPA(nlhdlrInitSepaSoc)
Definition: nlhdlr_soc.c:2581
SCIProwprepGetNVars
int SCIProwprepGetNVars(SCIP_ROWPREP *rowprep)
Definition: misc_rowprep.c:629
SCIP_DIDNOTFIND
Definition: type_result.h:44
SCIP_Cons
Definition: struct_cons.h:46
SCIPerrorMessage
#define SCIPerrorMessage
Definition: pub_message.h:64
SCIProwprepIsLocal
SCIP_Bool SCIProwprepIsLocal(SCIP_ROWPREP *rowprep)
Definition: misc_rowprep.c:679
SCIPisLT
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:458
buildQuadExprMatrix
static void buildQuadExprMatrix(SCIP *scip, SCIP_EXPR *quadexpr, SCIP_HASHMAP *expr2idx, int nexprs, SCIP_Real *quadmatrix, SCIP_Real *linvector)
Definition: nlhdlr_soc.c:990
SCIPnlhdlrGetName
const char * SCIPnlhdlrGetName(SCIP_NLHDLR *nlhdlr)
Definition: nlhdlr.c:166
SCIP_Conshdlr
Definition: struct_cons.h:126
SCIPgetVarExprVar
SCIP_VAR * SCIPgetVarExprVar(SCIP_EXPR *expr)
Definition: expr_var.c:416
SCIPhashsetFree
void SCIPhashsetFree(SCIP_HASHSET **hashset, BMS_BLKMEM *blkmem)
Definition: misc.c:3790
SCIPincludeNlhdlrNonlinear
SCIP_RETCODE SCIPincludeNlhdlrNonlinear(SCIP *scip, SCIP_NLHDLR **nlhdlr, const char *name, const char *desc, int detectpriority, int enfopriority, SCIP_DECL_NLHDLRDETECT((*detect)), SCIP_DECL_NLHDLREVALAUX((*evalaux)), SCIP_NLHDLRDATA *nlhdlrdata)
Definition: cons_nonlinear.c:15168
SCIPfreeBufferArrayNull
#define SCIPfreeBufferArrayNull(scip, ptr)
Definition: scip_mem.h:137
SCIPblkmem
BMS_BLKMEM * SCIPblkmem(SCIP *scip)
Definition: scip_mem.c:57
createDisaggrRow
static SCIP_RETCODE createDisaggrRow(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_EXPR *expr, SCIP_NLHDLREXPRDATA *nlhdlrexprdata)
Definition: nlhdlr_soc.c:282
SCIPvarGetName
const char * SCIPvarGetName(SCIP_VAR *var)
Definition: var.c:17420
nlhdlr_soc.h
soc nonlinear handler
SCIPhashmapFree
void SCIPhashmapFree(SCIP_HASHMAP **hashmap)
Definition: misc.c:3108
SCIPgetExprNonlinear
SCIP_EXPR * SCIPgetExprNonlinear(SCIP_CONS *cons)
Definition: cons_nonlinear.c:13705
expr_pow.h
power and signed power expression handlers
SCIP_OKAY
Definition: type_retcode.h:42
SCIPisExprSum
SCIP_Bool SCIPisExprSum(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1453
SCIPfreeRowprep
void SCIPfreeRowprep(SCIP *scip, SCIP_ROWPREP **rowprep)
Definition: misc_rowprep.c:583
SCIP_CALL
#define SCIP_CALL(x)
Definition: def.h:380
SCIProwprepAddSide
void SCIProwprepAddSide(SCIP_ROWPREP *rowprep, SCIP_Real side)
Definition: misc_rowprep.c:746
SCIP_DECL_NLHDLREVALAUX
static SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxSoc)
Definition: nlhdlr_soc.c:2487
SCIP_Interval::sup
SCIP_Real sup
Definition: intervalarith.h:56
SCIPverbMessage
void SCIPverbMessage(SCIP *scip, SCIP_VERBLEVEL msgverblevel, FILE *file, const char *formatstr,...)
Definition: scip_message.c:225
SCIPgetRowprepViolation
SCIP_Real SCIPgetRowprepViolation(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_SOL *sol, SCIP_Bool *reliable)
Definition: misc_rowprep.c:972
SCIPaddRow
SCIP_RETCODE SCIPaddRow(SCIP *scip, SCIP_ROW *row, SCIP_Bool forcecut, SCIP_Bool *infeasible)
Definition: scip_cut.c:250
SCIPnlhdlrSetSollinearize
void SCIPnlhdlrSetSollinearize(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRSOLLINEARIZE((*sollinearize)))
Definition: nlhdlr.c:154
SCIPdebugGetSolVal
#define SCIPdebugGetSolVal(scip, var, val)
Definition: debug.h:299
SCIPsortedvecFindInt
SCIP_Bool SCIPsortedvecFindInt(int *intarray, int val, int len, int *pos)
lapack_calls.h
interface methods for lapack functions
SCIPnlhdlrSetSepa
void SCIPnlhdlrSetSepa(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRINITSEPA((*initsepa)), SCIP_DECL_NLHDLRENFO((*enfo)), SCIP_DECL_NLHDLRESTIMATE((*estimate)), SCIP_DECL_NLHDLREXITSEPA((*exitsepa)))
Definition: nlhdlr.c:136
SCIP_Expr
Definition: struct_expr.h:105
SCIPgetExprNLocksNegNonlinear
int SCIPgetExprNLocksNegNonlinear(SCIP_EXPR *expr)
Definition: cons_nonlinear.c:14245
SCIPallocBufferArray
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:124
SCIP_Bool
#define SCIP_Bool
Definition: def.h:91
SCIP_DECL_NLHDLRDETECT
static SCIP_DECL_NLHDLRDETECT(nlhdlrDetectSoc)
Definition: nlhdlr_soc.c:2440
SCIPvarMarkRelaxationOnly
void SCIPvarMarkRelaxationOnly(SCIP_VAR *var)
Definition: var.c:17725
cons_nonlinear.h
constraint handler for nonlinear constraints specified by algebraic expressions
SCIP_DECL_NLHDLRCOPYHDLR
static SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrSoc)
Definition: nlhdlr_soc.c:2374
SCIP_DECL_NLHDLRINIT
#define SCIP_DECL_NLHDLRINIT(x)
Definition: type_nlhdlr.h:105
SCIP_Row
Definition: struct_lp.h:201
SCIPaddPoolCut
SCIP_RETCODE SCIPaddPoolCut(SCIP *scip, SCIP_ROW *row)
Definition: scip_cut.c:361
MIN
#define MIN(x, y)
Definition: def.h:243
debug.h
methods for debugging
SCIPgetCutEfficacy
SCIP_Real SCIPgetCutEfficacy(SCIP *scip, SCIP_SOL *sol, SCIP_ROW *cut)
Definition: scip_cut.c:94
nlhdlrExitSoc
#define nlhdlrExitSoc
Definition: nlhdlr_soc.c:2434
SCIPisExprPower
SCIP_Bool SCIPisExprPower(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1475
SCIPisCutApplicable
SCIP_Bool SCIPisCutApplicable(SCIP *scip, SCIP_ROW *cut)
Definition: scip_cut.c:207
SCIPisInfinity
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
Definition: scip_numerics.c:521
SCIP_LONGINT_FORMAT
#define SCIP_LONGINT_FORMAT
Definition: def.h:165
SCIPgetExprNLocksPosNonlinear
int SCIPgetExprNLocksPosNonlinear(SCIP_EXPR *expr)
Definition: cons_nonlinear.c:14234
generateCutSolSOC
static SCIP_RETCODE generateCutSolSOC(SCIP *scip, SCIP_ROWPREP **rowprep, SCIP_EXPR *expr, SCIP_CONS *cons, SCIP_NLHDLREXPRDATA *nlhdlrexprdata, SCIP_Real mincutviolation, SCIP_Real rhsval)
Definition: nlhdlr_soc.c:494
SCIPreleaseRow
SCIP_RETCODE SCIPreleaseRow(SCIP *scip, SCIP_ROW **row)
Definition: scip_lp.c:1562
MAX
#define MAX(x, y)
Definition: def.h:239
SCIPhashmapGetNElements
int SCIPhashmapGetNElements(SCIP_HASHMAP *hashmap)
Definition: misc.c:3533
SCIPisGT
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:484
SCIPisIntegral
SCIP_Bool SCIPisIntegral(SCIP *scip, SCIP_Real val)
Definition: scip_numerics.c:594
SCIPaddVar
SCIP_RETCODE SCIPaddVar(SCIP *scip, SCIP_VAR *var)
Definition: scip_prob.c:1668
evalSingleTerm
static SCIP_Real evalSingleTerm(SCIP *scip, SCIP_NLHDLREXPRDATA *nlhdlrexprdata, int k)
Definition: nlhdlr_soc.c:444
SCIPgetHugeValue
SCIP_Real SCIPgetHugeValue(SCIP *scip)
Definition: scip_numerics.c:547
NLHDLR_DESC
#define NLHDLR_DESC
Definition: nlhdlr_soc.c:57
NLHDLR_NAME
#define NLHDLR_NAME
Definition: nlhdlr_soc.c:56
SCIP_LOCKTYPE_MODEL
Definition: type_var.h:97
SCIPisExprVar
SCIP_Bool SCIPisExprVar(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1431
SCIPhashsetInsert
SCIP_RETCODE SCIPhashsetInsert(SCIP_HASHSET *hashset, BMS_BLKMEM *blkmem, void *element)
Definition: misc.c:3800
SCIPincludeNlhdlrSoc
SCIP_RETCODE SCIPincludeNlhdlrSoc(SCIP *scip)
Definition: nlhdlr_soc.c:3149
SCIPnlhdlrSetInitExit
void SCIPnlhdlrSetInitExit(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRINIT((*init)), SCIP_DECL_NLHDLREXIT((*exit_)))
Definition: nlhdlr.c:110
SCIP_Real
#define SCIP_Real
Definition: def.h:173
SCIP_DECL_NLHDLRENFO
static SCIP_DECL_NLHDLRENFO(nlhdlrEnfoSoc)
Definition: nlhdlr_soc.c:3001
SCIP_SIDETYPE_RIGHT
Definition: type_lp.h:65
detectSOC
static SCIP_RETCODE detectSOC(SCIP *scip, SCIP_NLHDLRDATA *nlhdlrdata, SCIP_EXPR *expr, SCIP_Real conslhs, SCIP_Real consrhs, SCIP_NLHDLREXPRDATA **nlhdlrexprdata, SCIP_Bool *enforcebelow, SCIP_Bool *success)
Definition: nlhdlr_soc.c:2324
SCIP_INVALID
#define SCIP_INVALID
Definition: def.h:193
SCIPnlhdlrGetData
SCIP_NLHDLRDATA * SCIPnlhdlrGetData(SCIP_NLHDLR *nlhdlr)
Definition: nlhdlr.c:216
SCIP_DECL_NLHDLREXIT
#define SCIP_DECL_NLHDLREXIT(x)
Definition: type_nlhdlr.h:114
SCIPgetRowSolFeasibility
SCIP_Real SCIPgetRowSolFeasibility(SCIP *scip, SCIP_ROW *row, SCIP_SOL *sol)
Definition: scip_lp.c:2167
SCIP_NLHDLREXPRDATA
struct SCIP_NlhdlrExprData SCIP_NLHDLREXPRDATA
Definition: type_nlhdlr.h:453
SCIPdebugAddSolVal
#define SCIPdebugAddSolVal(scip, var, val)
Definition: debug.h:298
SCIPcreateRowprep
SCIP_RETCODE SCIPcreateRowprep(SCIP *scip, SCIP_ROWPREP **rowprep, SCIP_SIDETYPE sidetype, SCIP_Bool local)
Definition: misc_rowprep.c:563
SCIP_NLHDLRDATA
struct SCIP_NlhdlrData SCIP_NLHDLRDATA
Definition: type_nlhdlr.h:452
SCIPisZero
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
Definition: scip_numerics.c:558
checkAndCollectQuadratic
static SCIP_RETCODE checkAndCollectQuadratic(SCIP *scip, SCIP_EXPR *quadexpr, SCIP_HASHMAP *expr2idx, SCIP_EXPR **occurringexprs, int *nexprs, SCIP_Bool *success)
Definition: nlhdlr_soc.c:875
expr_sum.h
sum expression handler
SCIPisLE
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:471
SCIPaddRowprepTerm
SCIP_RETCODE SCIPaddRowprepTerm(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_VAR *var, SCIP_Real coef)
Definition: misc_rowprep.c:913
SCIPfreeBlockMemoryArrayNull
#define SCIPfreeBlockMemoryArrayNull(scip, ptr, num)
Definition: scip_mem.h:111
SCIPisSOCNonlinear
SCIP_RETCODE SCIPisSOCNonlinear(SCIP *scip, SCIP_CONS *cons, SCIP_Bool compeigenvalues, SCIP_Bool *success, SCIP_SIDETYPE *sidetype, SCIP_VAR ***vars, SCIP_Real **offsets, SCIP_Real **transcoefs, int **transcoefsidx, int **termbegins, int *nvars, int *nterms)
Definition: nlhdlr_soc.c:3207
SCIPdismantleExpr
SCIP_RETCODE SCIPdismantleExpr(SCIP *scip, FILE *file, SCIP_EXPR *expr)
Definition: scip_expr.c:1608
SCIPgetExprAuxVarNonlinear
SCIP_VAR * SCIPgetExprAuxVarNonlinear(SCIP_EXPR *expr)
Definition: cons_nonlinear.c:14259
BMSclearMemoryArray
#define BMSclearMemoryArray(ptr, num)
Definition: memory.h:130
SCIPallocClearBlockMemory
#define SCIPallocClearBlockMemory(scip, ptr)
Definition: scip_mem.h:91
SCIP_DECL_NLHDLRFREEEXPRDATA
static SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataSoc)
Definition: nlhdlr_soc.c:2398
SCIPhashmapGetImageInt
int SCIPhashmapGetImageInt(SCIP_HASHMAP *hashmap, void *origin)
Definition: misc.c:3281
pub_nlhdlr.h
public functions of nonlinear handlers of nonlinear constraints
SCIPcreateEmptyRowConshdlr
SCIP_RETCODE SCIPcreateEmptyRowConshdlr(SCIP *scip, SCIP_ROW **row, SCIP_CONSHDLR *conshdlr, const char *name, SCIP_Real lhs, SCIP_Real rhs, SCIP_Bool local, SCIP_Bool modifiable, SCIP_Bool removable)
Definition: scip_lp.c:1391
SCIP_DECL_NLHDLREXITSEPA
static SCIP_DECL_NLHDLREXITSEPA(nlhdlrExitSepaSoc)
Definition: nlhdlr_soc.c:2983
SCIPgetNLPs
SCIP_Longint SCIPgetNLPs(SCIP *scip)
Definition: scip_solvingstats.c:427
SCIPABORT
#define SCIPABORT()
Definition: def.h:352
freeNlhdlrExprData
static SCIP_RETCODE freeNlhdlrExprData(SCIP *scip, SCIP_NLHDLREXPRDATA **nlhdlrexprdata)
Definition: nlhdlr_soc.c:384
SCIPround
SCIP_Real SCIPround(SCIP *scip, SCIP_Real val)
Definition: scip_numerics.c:655
SCIPgetLhsNonlinear
SCIP_Real SCIPgetLhsNonlinear(SCIP_CONS *cons)
Definition: cons_nonlinear.c:13721
SCIP_NLHDLR_METHOD_SEPABOTH
#define SCIP_NLHDLR_METHOD_SEPABOTH
Definition: type_nlhdlr.h:53
SCIPgetSolVal
SCIP_Real SCIPgetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var)
Definition: scip_sol.c:1217
SCIP_NLHDLR_METHOD_SEPAABOVE
#define SCIP_NLHDLR_METHOD_SEPAABOVE
Definition: type_nlhdlr.h:52
SCIPaddRealParam
SCIP_RETCODE SCIPaddRealParam(SCIP *scip, const char *name, const char *desc, SCIP_Real *valueptr, SCIP_Bool isadvanced, SCIP_Real defaultvalue, SCIP_Real minvalue, SCIP_Real maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:139
freeDisaggrVars
static SCIP_RETCODE freeDisaggrVars(SCIP *scip, SCIP_NLHDLREXPRDATA *nlhdlrexprdata)
Definition: nlhdlr_soc.c:251
SCIP_CUTOFF
Definition: type_result.h:48
SCIPaddBoolParam
SCIP_RETCODE SCIPaddBoolParam(SCIP *scip, const char *name, const char *desc, SCIP_Bool *valueptr, SCIP_Bool isadvanced, SCIP_Bool defaultvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:57
SCIPlapackIsAvailable
SCIP_Bool SCIPlapackIsAvailable(void)
Definition: lapack_calls.c:121
SCIP_RowPrep
Definition: struct_misc.h:286
SCIPnlhdlrSetCopyHdlr
void SCIPnlhdlrSetCopyHdlr(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRCOPYHDLR((*copy)))
Definition: nlhdlr.c:76
SCIP_SIDETYPE
enum SCIP_SideType SCIP_SIDETYPE
Definition: type_lp.h:67
SCIP_VARTYPE_CONTINUOUS
Definition: type_var.h:71