Scippy

SCIP

Solving Constraint Integer Programs

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