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