1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file scip_nonlinear.c
17 * @ingroup OTHER_CFILES
18 * @brief public methods for nonlinear functions
19 * @author Tobias Achterberg
20 * @author Timo Berthold
21 * @author Gerald Gamrath
22 * @author Leona Gottwald
23 * @author Stefan Heinz
24 * @author Gregor Hendel
25 * @author Thorsten Koch
26 * @author Alexander Martin
27 * @author Marc Pfetsch
28 * @author Michael Winkler
29 * @author Kati Wolter
30 *
31 * @todo check all SCIP_STAGE_* switches, and include the new stages TRANSFORMED and INITSOLVE
32 */
33
34 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
35
36 #include "blockmemshell/memory.h"
37 #include "nlpi/nlpi.h"
38 #include "nlpi/pub_expr.h"
39 #include "nlpi/type_expr.h"
40 #include "scip/dbldblarith.h"
41 #include "scip/pub_lp.h"
42 #include "scip/pub_message.h"
43 #include "scip/pub_misc.h"
44 #include "scip/pub_nlp.h"
45 #include "scip/pub_var.h"
46 #include "scip/scip_mem.h"
47 #include "scip/scip_message.h"
48 #include "scip/scip_nonlinear.h"
49 #include "scip/scip_numerics.h"
50 #include "scip/scip_prob.h"
51
52 /** computes coefficients of linearization of a square term in a reference point */
SCIPaddSquareLinearization(SCIP * scip,SCIP_Real sqrcoef,SCIP_Real refpoint,SCIP_Bool isint,SCIP_Real * lincoef,SCIP_Real * linconstant,SCIP_Bool * success)53 void SCIPaddSquareLinearization(
54 SCIP* scip, /**< SCIP data structure */
55 SCIP_Real sqrcoef, /**< coefficient of square term */
56 SCIP_Real refpoint, /**< point where to linearize */
57 SCIP_Bool isint, /**< whether corresponding variable is a discrete variable, and thus linearization could be moved */
58 SCIP_Real* lincoef, /**< buffer to add coefficient of linearization */
59 SCIP_Real* linconstant, /**< buffer to add constant of linearization */
60 SCIP_Bool* success /**< buffer to set to FALSE if linearization has failed due to large numbers */
61 )
62 {
63 assert(scip != NULL);
64 assert(lincoef != NULL);
65 assert(linconstant != NULL);
66 assert(success != NULL);
67
68 if( sqrcoef == 0.0 )
69 return;
70
71 if( SCIPisInfinity(scip, REALABS(refpoint)) )
72 {
73 *success = FALSE;
74 return;
75 }
76
77 if( !isint || SCIPisIntegral(scip, refpoint) )
78 {
79 SCIP_Real tmp;
80
81 /* sqrcoef * x^2 -> tangent in refpoint = sqrcoef * 2 * refpoint * (x - refpoint) */
82
83 tmp = sqrcoef * refpoint;
84
85 if( SCIPisInfinity(scip, 2.0 * REALABS(tmp)) )
86 {
87 *success = FALSE;
88 return;
89 }
90
91 *lincoef += 2.0 * tmp;
92 tmp *= refpoint;
93 *linconstant -= tmp;
94 }
95 else
96 {
97 /* sqrcoef * x^2 -> secant between f=floor(refpoint) and f+1 = sqrcoef * (f^2 + ((f+1)^2 - f^2) * (x-f))
98 * = sqrcoef * (-f*(f+1) + (2*f+1)*x)
99 */
100 SCIP_Real f;
101 SCIP_Real coef;
102 SCIP_Real constant;
103
104 f = SCIPfloor(scip, refpoint);
105
106 coef = sqrcoef * (2.0 * f + 1.0);
107 constant = -sqrcoef * f * (f + 1.0);
108
109 if( SCIPisInfinity(scip, REALABS(coef)) || SCIPisInfinity(scip, REALABS(constant)) )
110 {
111 *success = FALSE;
112 return;
113 }
114
115 *lincoef += coef;
116 *linconstant += constant;
117 }
118 }
119
120 /** computes coefficients of secant of a square term */
SCIPaddSquareSecant(SCIP * scip,SCIP_Real sqrcoef,SCIP_Real lb,SCIP_Real ub,SCIP_Real refpoint,SCIP_Real * lincoef,SCIP_Real * linconstant,SCIP_Bool * success)121 void SCIPaddSquareSecant(
122 SCIP* scip, /**< SCIP data structure */
123 SCIP_Real sqrcoef, /**< coefficient of square term */
124 SCIP_Real lb, /**< lower bound on variable */
125 SCIP_Real ub, /**< upper bound on variable */
126 SCIP_Real refpoint, /**< point for which to compute value of linearization */
127 SCIP_Real* lincoef, /**< buffer to add coefficient of secant */
128 SCIP_Real* linconstant, /**< buffer to add constant of secant */
129 SCIP_Bool* success /**< buffer to set to FALSE if secant has failed due to large numbers or unboundedness */
130 )
131 {
132 SCIP_Real coef;
133 SCIP_Real constant;
134
135 assert(scip != NULL);
136 assert(!SCIPisInfinity(scip, lb));
137 assert(!SCIPisInfinity(scip, -ub));
138 assert(SCIPisLE(scip, lb, ub));
139 assert(SCIPisLE(scip, lb, refpoint));
140 assert(SCIPisGE(scip, ub, refpoint));
141 assert(lincoef != NULL);
142 assert(linconstant != NULL);
143 assert(success != NULL);
144
145 if( sqrcoef == 0.0 )
146 return;
147
148 if( SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub) )
149 {
150 /* unboundedness */
151 *success = FALSE;
152 return;
153 }
154
155 /* sqrcoef * x^2 -> sqrcoef * (lb * lb + (ub*ub - lb*lb)/(ub-lb) * (x-lb)) = sqrcoef * (lb*lb + (ub+lb)*(x-lb))
156 * = sqrcoef * ((lb+ub)*x - lb*ub)
157 */
158 coef = sqrcoef * (lb + ub);
159 constant = -sqrcoef * lb * ub;
160 if( SCIPisInfinity(scip, REALABS(coef)) || SCIPisInfinity(scip, REALABS(constant)) )
161 {
162 *success = FALSE;
163 return;
164 }
165
166 *lincoef += coef;
167 *linconstant += constant;
168 }
169
170 /** computes coefficients of linearization of a bilinear term in a reference point */
SCIPaddBilinLinearization(SCIP * scip,SCIP_Real bilincoef,SCIP_Real refpointx,SCIP_Real refpointy,SCIP_Real * lincoefx,SCIP_Real * lincoefy,SCIP_Real * linconstant,SCIP_Bool * success)171 void SCIPaddBilinLinearization(
172 SCIP* scip, /**< SCIP data structure */
173 SCIP_Real bilincoef, /**< coefficient of bilinear term */
174 SCIP_Real refpointx, /**< point where to linearize first variable */
175 SCIP_Real refpointy, /**< point where to linearize second variable */
176 SCIP_Real* lincoefx, /**< buffer to add coefficient of first variable in linearization */
177 SCIP_Real* lincoefy, /**< buffer to add coefficient of second variable in linearization */
178 SCIP_Real* linconstant, /**< buffer to add constant of linearization */
179 SCIP_Bool* success /**< buffer to set to FALSE if linearization has failed due to large numbers */
180 )
181 {
182 SCIP_Real constant;
183
184 assert(scip != NULL);
185 assert(lincoefx != NULL);
186 assert(lincoefy != NULL);
187 assert(linconstant != NULL);
188 assert(success != NULL);
189
190 if( bilincoef == 0.0 )
191 return;
192
193 if( SCIPisInfinity(scip, REALABS(refpointx)) || SCIPisInfinity(scip, REALABS(refpointy)) )
194 {
195 *success = FALSE;
196 return;
197 }
198
199 /* bilincoef * x * y -> bilincoef * (refpointx * refpointy + refpointy * (x - refpointx) + refpointx * (y - refpointy))
200 * = -bilincoef * refpointx * refpointy + bilincoef * refpointy * x + bilincoef * refpointx * y
201 */
202
203 constant = -bilincoef * refpointx * refpointy;
204
205 if( SCIPisInfinity(scip, REALABS(bilincoef * refpointx)) || SCIPisInfinity(scip, REALABS(bilincoef * refpointy))
206 || SCIPisInfinity(scip, REALABS(constant)) )
207 {
208 *success = FALSE;
209 return;
210 }
211
212 *lincoefx += bilincoef * refpointy;
213 *lincoefy += bilincoef * refpointx;
214 *linconstant += constant;
215 }
216
217 /** computes coefficients of McCormick under- or overestimation of a bilinear term */
SCIPaddBilinMcCormick(SCIP * scip,SCIP_Real bilincoef,SCIP_Real lbx,SCIP_Real ubx,SCIP_Real refpointx,SCIP_Real lby,SCIP_Real uby,SCIP_Real refpointy,SCIP_Bool overestimate,SCIP_Real * lincoefx,SCIP_Real * lincoefy,SCIP_Real * linconstant,SCIP_Bool * success)218 void SCIPaddBilinMcCormick(
219 SCIP* scip, /**< SCIP data structure */
220 SCIP_Real bilincoef, /**< coefficient of bilinear term */
221 SCIP_Real lbx, /**< lower bound on first variable */
222 SCIP_Real ubx, /**< upper bound on first variable */
223 SCIP_Real refpointx, /**< reference point for first variable */
224 SCIP_Real lby, /**< lower bound on second variable */
225 SCIP_Real uby, /**< upper bound on second variable */
226 SCIP_Real refpointy, /**< reference point for second variable */
227 SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
228 SCIP_Real* lincoefx, /**< buffer to add coefficient of first variable in linearization */
229 SCIP_Real* lincoefy, /**< buffer to add coefficient of second variable in linearization */
230 SCIP_Real* linconstant, /**< buffer to add constant of linearization */
231 SCIP_Bool* success /**< buffer to set to FALSE if linearization has failed due to large numbers */
232 )
233 {
234 SCIP_Real constant;
235 SCIP_Real coefx;
236 SCIP_Real coefy;
237
238 assert(scip != NULL);
239 assert(!SCIPisInfinity(scip, lbx));
240 assert(!SCIPisInfinity(scip, -ubx));
241 assert(!SCIPisInfinity(scip, lby));
242 assert(!SCIPisInfinity(scip, -uby));
243 assert(SCIPisLE(scip, lbx, ubx));
244 assert(SCIPisLE(scip, lby, uby));
245 assert(SCIPisLE(scip, lbx, refpointx));
246 assert(SCIPisGE(scip, ubx, refpointx));
247 assert(SCIPisLE(scip, lby, refpointy));
248 assert(SCIPisGE(scip, uby, refpointy));
249 assert(lincoefx != NULL);
250 assert(lincoefy != NULL);
251 assert(linconstant != NULL);
252 assert(success != NULL);
253
254 if( bilincoef == 0.0 )
255 return;
256
257 if( overestimate )
258 bilincoef = -bilincoef;
259
260 if( SCIPisRelEQ(scip, lbx, ubx) && SCIPisRelEQ(scip, lby, uby) )
261 {
262 /* both x and y are mostly fixed */
263 SCIP_Real cand1;
264 SCIP_Real cand2;
265 SCIP_Real cand3;
266 SCIP_Real cand4;
267
268 coefx = 0.0;
269 coefy = 0.0;
270
271 /* estimate x * y by constant */
272 cand1 = lbx * lby;
273 cand2 = lbx * uby;
274 cand3 = ubx * lby;
275 cand4 = ubx * uby;
276
277 /* take most conservative value for underestimator */
278 if( bilincoef < 0.0 )
279 constant = bilincoef * MAX( MAX(cand1, cand2), MAX(cand3, cand4) );
280 else
281 constant = bilincoef * MIN( MIN(cand1, cand2), MIN(cand3, cand4) );
282 }
283 else if( bilincoef > 0.0 )
284 {
285 /* either x or y is not fixed and coef > 0.0 */
286 if( !SCIPisInfinity(scip, -lbx) && !SCIPisInfinity(scip, -lby) &&
287 (SCIPisInfinity(scip, ubx) || SCIPisInfinity(scip, uby)
288 || (uby - refpointy) * (ubx - refpointx) >= (refpointy - lby) * (refpointx - lbx)) )
289 {
290 if( SCIPisRelEQ(scip, lbx, ubx) )
291 {
292 /* x*y = lbx * y + (x-lbx) * y >= lbx * y + (x-lbx) * lby >= lbx * y + min{(ubx-lbx) * lby, 0 * lby} */
293 coefx = 0.0;
294 coefy = bilincoef * lbx;
295 constant = bilincoef * (lby < 0.0 ? (ubx-lbx) * lby : 0.0);
296 }
297 else if( SCIPisRelEQ(scip, lby, uby) )
298 {
299 /* x*y = lby * x + (y-lby) * x >= lby * x + (y-lby) * lbx >= lby * x + min{(uby-lby) * lbx, 0 * lbx} */
300 coefx = bilincoef * lby;
301 coefy = 0.0;
302 constant = bilincoef * (lbx < 0.0 ? (uby-lby) * lbx : 0.0);
303 }
304 else
305 {
306 coefx = bilincoef * lby;
307 coefy = bilincoef * lbx;
308 constant = -bilincoef * lbx * lby;
309 }
310 }
311 else if( !SCIPisInfinity(scip, ubx) && !SCIPisInfinity(scip, uby) )
312 {
313 if( SCIPisRelEQ(scip, lbx, ubx) )
314 {
315 /* x*y = ubx * y + (x-ubx) * y >= ubx * y + (x-ubx) * uby >= ubx * y + min{(lbx-ubx) * uby, 0 * uby} */
316 coefx = 0.0;
317 coefy = bilincoef * ubx;
318 constant = bilincoef * (uby > 0.0 ? (lbx - ubx) * uby : 0.0);
319 }
320 else if( SCIPisRelEQ(scip, lby, uby) )
321 {
322 /* x*y = uby * x + (y-uby) * x >= uby * x + (y-uby) * ubx >= uby * x + min{(lby-uby) * ubx, 0 * ubx} */
323 coefx = bilincoef * uby;
324 coefy = 0.0;
325 constant = bilincoef * (ubx > 0.0 ? (lby - uby) * ubx : 0.0);
326 }
327 else
328 {
329 coefx = bilincoef * uby;
330 coefy = bilincoef * ubx;
331 constant = -bilincoef * ubx * uby;
332 }
333 }
334 else
335 {
336 *success = FALSE;
337 return;
338 }
339 }
340 else
341 {
342 /* either x or y is not fixed and coef < 0.0 */
343 if( !SCIPisInfinity(scip, ubx) && !SCIPisInfinity(scip, -lby) &&
344 (SCIPisInfinity(scip, -lbx) || SCIPisInfinity(scip, uby)
345 || (ubx - lbx) * (refpointy - lby) <= (uby - lby) * (refpointx - lbx)) )
346 {
347 if( SCIPisRelEQ(scip, lbx, ubx) )
348 {
349 /* x*y = ubx * y + (x-ubx) * y <= ubx * y + (x-ubx) * lby <= ubx * y + max{(lbx-ubx) * lby, 0 * lby} */
350 coefx = 0.0;
351 coefy = bilincoef * ubx;
352 constant = bilincoef * (lby < 0.0 ? (lbx - ubx) * lby : 0.0);
353 }
354 else if( SCIPisRelEQ(scip, lby, uby) )
355 {
356 /* x*y = lby * x + (y-lby) * x <= lby * x + (y-lby) * ubx <= lby * x + max{(uby-lby) * ubx, 0 * ubx} */
357 coefx = bilincoef * lby;
358 coefy = 0.0;
359 constant = bilincoef * (ubx > 0.0 ? (uby - lby) * ubx : 0.0);
360 }
361 else
362 {
363 coefx = bilincoef * lby;
364 coefy = bilincoef * ubx;
365 constant = -bilincoef * ubx * lby;
366 }
367 }
368 else if( !SCIPisInfinity(scip, -lbx) && !SCIPisInfinity(scip, uby) )
369 {
370 if( SCIPisRelEQ(scip, lbx, ubx) )
371 {
372 /* x*y = lbx * y + (x-lbx) * y <= lbx * y + (x-lbx) * uby <= lbx * y + max{(ubx-lbx) * uby, 0 * uby} */
373 coefx = 0.0;
374 coefy = bilincoef * lbx;
375 constant = bilincoef * (uby > 0.0 ? (ubx - lbx) * uby : 0.0);
376 }
377 else if( SCIPisRelEQ(scip, lby, uby) )
378 {
379 /* x*y = uby * x + (y-uby) * x <= uby * x + (y-uby) * lbx <= uby * x + max{(lby-uby) * lbx, 0 * lbx} */
380 coefx = bilincoef * uby;
381 coefy = 0.0;
382 constant = bilincoef * (lbx < 0.0 ? (lby - uby) * lbx : 0.0);
383 }
384 else
385 {
386 coefx = bilincoef * uby;
387 coefy = bilincoef * lbx;
388 constant = -bilincoef * lbx * uby;
389 }
390 }
391 else
392 {
393 *success = FALSE;
394 return;
395 }
396 }
397
398 if( SCIPisInfinity(scip, REALABS(coefx)) || SCIPisInfinity(scip, REALABS(coefy))
399 || SCIPisInfinity(scip, REALABS(constant)) )
400 {
401 *success = FALSE;
402 return;
403 }
404
405 if( overestimate )
406 {
407 coefx = -coefx;
408 coefy = -coefy;
409 constant = -constant;
410 }
411
412 SCIPdebugMsg(scip, "%.15g * x[%.15g,%.15g] * y[%.15g,%.15g] %c= %.15g * x %+.15g * y %+.15g\n", bilincoef, lbx, ubx,
413 lby, uby, overestimate ? '<' : '>', coefx, coefy, constant);
414
415 *lincoefx += coefx;
416 *lincoefy += coefy;
417 *linconstant += constant;
418 }
419
420
421 /** computes coefficients of linearization of a bilinear term in a reference point when given a linear inequality
422 * involving only the variables of the bilinear term
423 *
424 * @note the formulas are extracted from "Convex envelopes of bivariate functions through the solution of KKT systems"
425 * by Marco Locatelli
426 */
SCIPcomputeBilinEnvelope1(SCIP * scip,SCIP_Real bilincoef,SCIP_Real lbx,SCIP_Real ubx,SCIP_Real refpointx,SCIP_Real lby,SCIP_Real uby,SCIP_Real refpointy,SCIP_Bool overestimate,SCIP_Real xcoef,SCIP_Real ycoef,SCIP_Real constant,SCIP_Real * RESTRICT lincoefx,SCIP_Real * RESTRICT lincoefy,SCIP_Real * RESTRICT linconstant,SCIP_Bool * RESTRICT success)427 void SCIPcomputeBilinEnvelope1(
428 SCIP* scip, /**< SCIP data structure */
429 SCIP_Real bilincoef, /**< coefficient of bilinear term */
430 SCIP_Real lbx, /**< lower bound on first variable */
431 SCIP_Real ubx, /**< upper bound on first variable */
432 SCIP_Real refpointx, /**< reference point for first variable */
433 SCIP_Real lby, /**< lower bound on second variable */
434 SCIP_Real uby, /**< upper bound on second variable */
435 SCIP_Real refpointy, /**< reference point for second variable */
436 SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
437 SCIP_Real xcoef, /**< x coefficient of linear inequality; must be in {-1,0,1} */
438 SCIP_Real ycoef, /**< y coefficient of linear inequality */
439 SCIP_Real constant, /**< constant of linear inequality */
440 SCIP_Real* RESTRICT lincoefx, /**< buffer to store coefficient of first variable in linearization */
441 SCIP_Real* RESTRICT lincoefy, /**< buffer to store coefficient of second variable in linearization */
442 SCIP_Real* RESTRICT linconstant, /**< buffer to store constant of linearization */
443 SCIP_Bool* RESTRICT success /**< buffer to store whether linearization was successful */
444 )
445 {
446 SCIP_Real xs[2] = {lbx, ubx};
447 SCIP_Real ys[2] = {lby, uby};
448 SCIP_Real minx;
449 SCIP_Real maxx;
450 SCIP_Real miny;
451 SCIP_Real maxy;
452 SCIP_Real QUAD(lincoefyq);
453 SCIP_Real QUAD(lincoefxq);
454 SCIP_Real QUAD(linconstantq);
455 SCIP_Real QUAD(denomq);
456 SCIP_Real QUAD(mjq);
457 SCIP_Real QUAD(qjq);
458 SCIP_Real QUAD(xjq);
459 SCIP_Real QUAD(yjq);
460 SCIP_Real QUAD(tmpq);
461 SCIP_Real vx;
462 SCIP_Real vy;
463 int n;
464 int i;
465
466 assert(scip != NULL);
467 assert(!SCIPisInfinity(scip, lbx));
468 assert(!SCIPisInfinity(scip, -ubx));
469 assert(!SCIPisInfinity(scip, lby));
470 assert(!SCIPisInfinity(scip, -uby));
471 assert(SCIPisLE(scip, lbx, ubx));
472 assert(SCIPisLE(scip, lby, uby));
473 assert(SCIPisLE(scip, lbx, refpointx));
474 assert(SCIPisGE(scip, ubx, refpointx));
475 assert(SCIPisLE(scip, lby, refpointy));
476 assert(SCIPisGE(scip, uby, refpointy));
477 assert(lincoefx != NULL);
478 assert(lincoefy != NULL);
479 assert(linconstant != NULL);
480 assert(success != NULL);
481 assert(xcoef == 0.0 || xcoef == -1.0 || xcoef == 1.0); /*lint !e777*/
482 assert(ycoef != SCIP_INVALID && ycoef != 0.0); /*lint !e777*/
483 assert(constant != SCIP_INVALID); /*lint !e777*/
484
485 *success = FALSE;
486 *lincoefx = SCIP_INVALID;
487 *lincoefy = SCIP_INVALID;
488 *linconstant = SCIP_INVALID;
489
490 /* reference point does not satisfy linear inequality */
491 if( SCIPisFeasGT(scip, xcoef * refpointx - ycoef * refpointy - constant, 0.0) )
492 return;
493
494 /* compute minimal and maximal bounds on x and y for accepting the reference point */
495 minx = lbx + 0.01 * (ubx-lbx);
496 maxx = ubx - 0.01 * (ubx-lbx);
497 miny = lby + 0.01 * (uby-lby);
498 maxy = uby - 0.01 * (uby-lby);
499
500 /* check whether the reference point is in [minx,maxx]x[miny,maxy] */
501 if( SCIPisLE(scip, refpointx, minx) || SCIPisGE(scip, refpointx, maxx)
502 || SCIPisLE(scip, refpointy, miny) || SCIPisGE(scip, refpointy, maxy) )
503 return;
504
505 /* always consider xy without the bilinear coefficient */
506 if( bilincoef < 0.0 )
507 overestimate = !overestimate;
508
509 /* we use same notation as in "Convex envelopes of bivariate functions through the solution of KKT systems", 2016 */
510 /* mj = xcoef / ycoef */
511 SCIPquadprecDivDD(mjq, xcoef, ycoef);
512
513 /* qj = -constant / ycoef */
514 SCIPquadprecDivDD(qjq, -constant, ycoef);
515
516 /* mj > 0 => underestimate; mj < 0 => overestimate */
517 if( SCIPisNegative(scip, QUAD_TO_DBL(mjq)) != overestimate )
518 return;
519
520 /* get the corner point that satisfies the linear inequality xcoef*x <= ycoef*y + constant */
521 if( !overestimate )
522 {
523 ys[0] = uby;
524 ys[1] = lby;
525 }
526
527 vx = SCIP_INVALID;
528 vy = SCIP_INVALID;
529 n = 0;
530 for( i = 0; i < 2; ++i )
531 {
532 SCIP_Real activity = xcoef * xs[i] - ycoef * ys[i] - constant;
533 if( SCIPisLE(scip, activity, 0.0) )
534 {
535 /* corner point is satisfies inequality */
536 vx = xs[i];
537 vy = ys[i];
538 }
539 else if( SCIPisFeasGT(scip, activity, 0.0) )
540 /* corner point is clearly cut off */
541 ++n;
542 }
543
544 /* skip if no corner point satisfies the inequality or if no corner point is cut off (that is, all corner points satisfy the inequality almost [1e-9..1e-6]) */
545 if( n != 1 || vx == SCIP_INVALID || vy == SCIP_INVALID ) /*lint !e777*/
546 return;
547
548 /* denom = mj*(refpointx - vx) + vy - refpointy */
549 SCIPquadprecSumDD(denomq, refpointx, -vx); /* refpoint - vx */
550 SCIPquadprecProdQQ(denomq, denomq, mjq); /* mj * (refpoint - vx) */
551 SCIPquadprecSumQD(denomq, denomq, vy); /* mj * (refpoint - vx) + vy */
552 SCIPquadprecSumQD(denomq, denomq, -refpointy); /* mj * (refpoint - vx) + vy - refpointy */
553
554 if( SCIPisZero(scip, QUAD_TO_DBL(denomq)) )
555 return;
556
557 /* (xj,yj) is the projection onto the line xcoef*x = ycoef*y + constant */
558 /* xj = (refpointx*(vy - qj) - vx*(refpointy - qj)) / denom */
559 SCIPquadprecProdQD(xjq, qjq, -1.0); /* - qj */
560 SCIPquadprecSumQD(xjq, xjq, vy); /* vy - qj */
561 SCIPquadprecProdQD(xjq, xjq, refpointx); /* refpointx * (vy - qj) */
562 SCIPquadprecProdQD(tmpq, qjq, -1.0); /* - qj */
563 SCIPquadprecSumQD(tmpq, tmpq, refpointy); /* refpointy - qj */
564 SCIPquadprecProdQD(tmpq, tmpq, -vx); /* - vx * (refpointy - qj) */
565 SCIPquadprecSumQQ(xjq, xjq, tmpq); /* refpointx * (vy - qj) - vx * (refpointy - qj) */
566 SCIPquadprecDivQQ(xjq, xjq, denomq); /* (refpointx * (vy - qj) - vx * (refpointy - qj)) / denom */
567
568 /* yj = mj * xj + qj */
569 SCIPquadprecProdQQ(yjq, mjq, xjq);
570 SCIPquadprecSumQQ(yjq, yjq, qjq);
571
572 assert(SCIPisFeasEQ(scip, xcoef*QUAD_TO_DBL(xjq) - ycoef*QUAD_TO_DBL(yjq) - constant, 0.0));
573
574 /* check whether the projection is in [minx,maxx] x [miny,maxy]; this avoids numerical difficulties when the
575 * projection is close to the variable bounds
576 */
577 if( SCIPisLE(scip, QUAD_TO_DBL(xjq), minx) || SCIPisGE(scip, QUAD_TO_DBL(xjq), maxx)
578 || SCIPisLE(scip, QUAD_TO_DBL(yjq), miny) || SCIPisGE(scip, QUAD_TO_DBL(yjq), maxy) )
579 return;
580
581 assert(vy - QUAD_TO_DBL(mjq)*vx - QUAD_TO_DBL(qjq) != 0.0);
582
583 /* lincoefy = (mj*SQR(xj) - 2.0*mj*vx*xj - qj*vx + vx*vy) / (vy - mj*vx - qj) */
584 SCIPquadprecSquareQ(lincoefyq, xjq); /* xj^2 */
585 SCIPquadprecProdQQ(lincoefyq, lincoefyq, mjq); /* mj * xj^2 */
586 SCIPquadprecProdQQ(tmpq, mjq, xjq); /* mj * xj */
587 SCIPquadprecProdQD(tmpq, tmpq, -2.0 * vx); /* -2 * vx * mj * xj */
588 SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj */
589 SCIPquadprecProdQD(tmpq, qjq, -vx); /* -qj * vx */
590 SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj -qj * vx */
591 SCIPquadprecProdDD(tmpq, vx, vy); /* vx * vy */
592 SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj -qj * vx + vx * vy */
593 SCIPquadprecProdQD(tmpq, mjq, vx); /* mj * vx */
594 SCIPquadprecSumQD(tmpq, tmpq, -vy); /* -vy + mj * vx */
595 SCIPquadprecSumQQ(tmpq, tmpq, qjq); /* -vy + mj * vx + qj */
596 QUAD_SCALE(tmpq, -1.0); /* vy - mj * vx - qj */
597 SCIPquadprecDivQQ(lincoefyq, lincoefyq, tmpq); /* (mj * xj^2 -2 * vx * mj * xj -qj * vx + vx * vy) / (vy - mj * vx - qj) */
598
599 /* lincoefx = 2.0*mj*xj + qj - mj*(*lincoefy) */
600 SCIPquadprecProdQQ(lincoefxq, mjq, xjq); /* mj * xj */
601 QUAD_SCALE(lincoefxq, 2.0); /* 2 * mj * xj */
602 SCIPquadprecSumQQ(lincoefxq, lincoefxq, qjq); /* 2 * mj * xj + qj */
603 SCIPquadprecProdQQ(tmpq, mjq, lincoefyq); /* mj * lincoefy */
604 QUAD_SCALE(tmpq, -1.0); /* - mj * lincoefy */
605 SCIPquadprecSumQQ(lincoefxq, lincoefxq, tmpq); /* 2 * mj * xj + qj - mj * lincoefy */
606
607 /* linconstant = -mj*SQR(xj) - (*lincoefy)*qj */
608 SCIPquadprecSquareQ(linconstantq, xjq); /* xj^2 */
609 SCIPquadprecProdQQ(linconstantq, linconstantq, mjq); /* mj * xj^2 */
610 QUAD_SCALE(linconstantq, -1.0); /* - mj * xj^2 */
611 SCIPquadprecProdQQ(tmpq, lincoefyq, qjq); /* lincoefy * qj */
612 QUAD_SCALE(tmpq, -1.0); /* - lincoefy * qj */
613 SCIPquadprecSumQQ(linconstantq, linconstantq, tmpq); /* - mj * xj^2 - lincoefy * qj */
614
615 /* consider the bilinear coefficient */
616 SCIPquadprecProdQD(lincoefxq, lincoefxq, bilincoef);
617 SCIPquadprecProdQD(lincoefyq, lincoefyq, bilincoef);
618 SCIPquadprecProdQD(linconstantq, linconstantq, bilincoef);
619 *lincoefx = QUAD_TO_DBL(lincoefxq);
620 *lincoefy = QUAD_TO_DBL(lincoefyq);
621 *linconstant = QUAD_TO_DBL(linconstantq);
622
623 /* cut needs to be tight at (vx,vy) and (xj,yj); otherwise we consider the cut to be numerically bad */
624 *success = SCIPisFeasEQ(scip, (*lincoefx)*vx + (*lincoefy)*vy + (*linconstant), bilincoef*vx*vy)
625 && SCIPisFeasEQ(scip, (*lincoefx)*QUAD_TO_DBL(xjq) + (*lincoefy)*QUAD_TO_DBL(yjq) + (*linconstant), bilincoef*QUAD_TO_DBL(xjq)*QUAD_TO_DBL(yjq));
626
627 #ifndef NDEBUG
628 {
629 SCIP_Real activity = (*lincoefx)*refpointx + (*lincoefy)*refpointy + (*linconstant);
630
631 /* cut needs to under- or overestimate the bilinear term at the reference point */
632 if( bilincoef < 0.0 )
633 overestimate = !overestimate;
634
635 if( overestimate )
636 assert(SCIPisFeasGE(scip, activity, bilincoef*refpointx*refpointy));
637 else
638 assert(SCIPisFeasLE(scip, activity, bilincoef*refpointx*refpointy));
639 }
640 #endif
641 }
642
643 /** helper function to compute the convex envelope of a bilinear term when two linear inequalities are given; we
644 * use the same notation and formulas as in Locatelli 2016
645 */
646 static
computeBilinEnvelope2(SCIP * scip,SCIP_Real x,SCIP_Real y,SCIP_Real mi,SCIP_Real qi,SCIP_Real mj,SCIP_Real qj,SCIP_Real * RESTRICT xi,SCIP_Real * RESTRICT yi,SCIP_Real * RESTRICT xj,SCIP_Real * RESTRICT yj,SCIP_Real * RESTRICT xcoef,SCIP_Real * RESTRICT ycoef,SCIP_Real * RESTRICT constant)647 void computeBilinEnvelope2(
648 SCIP* scip, /**< SCIP data structure */
649 SCIP_Real x, /**< reference point for x */
650 SCIP_Real y, /**< reference point for y */
651 SCIP_Real mi, /**< coefficient of x in the first linear inequality */
652 SCIP_Real qi, /**< constant in the first linear inequality */
653 SCIP_Real mj, /**< coefficient of x in the second linear inequality */
654 SCIP_Real qj, /**< constant in the second linear inequality */
655 SCIP_Real* RESTRICT xi, /**< buffer to store x coordinate of the first point */
656 SCIP_Real* RESTRICT yi, /**< buffer to store y coordinate of the first point */
657 SCIP_Real* RESTRICT xj, /**< buffer to store x coordinate of the second point */
658 SCIP_Real* RESTRICT yj, /**< buffer to store y coordinate of the second point */
659 SCIP_Real* RESTRICT xcoef, /**< buffer to store the x coefficient of the envelope */
660 SCIP_Real* RESTRICT ycoef, /**< buffer to store the y coefficient of the envelope */
661 SCIP_Real* RESTRICT constant /**< buffer to store the constant of the envelope */
662 )
663 {
664 SCIP_Real QUAD(xiq);
665 SCIP_Real QUAD(yiq);
666 SCIP_Real QUAD(xjq);
667 SCIP_Real QUAD(yjq);
668 SCIP_Real QUAD(xcoefq);
669 SCIP_Real QUAD(ycoefq);
670 SCIP_Real QUAD(constantq);
671 SCIP_Real QUAD(tmpq);
672
673 assert(xi != NULL);
674 assert(yi != NULL);
675 assert(xj != NULL);
676 assert(yj != NULL);
677 assert(xcoef != NULL);
678 assert(ycoef != NULL);
679 assert(constant != NULL);
680
681 if( SCIPisEQ(scip, mi, mj) )
682 {
683 /* xi = (x + mi * y - qi) / (2.0*mi) */
684 SCIPquadprecProdDD(xiq, mi, y);
685 SCIPquadprecSumQD(xiq, xiq, x);
686 SCIPquadprecSumQD(xiq, xiq, -qi);
687 SCIPquadprecDivQD(xiq, xiq, 2.0 * mi);
688 assert(EPSEQ((x + mi * y - qi) / (2.0*mi), QUAD_TO_DBL(xiq), 1e-3));
689
690 /* yi = mi*(*xi) + qi */
691 SCIPquadprecProdQD(yiq, xiq, mi);
692 SCIPquadprecSumQD(yiq, yiq, qi);
693 assert(EPSEQ(mi*QUAD_TO_DBL(xiq) + qi, QUAD_TO_DBL(yiq), 1e-3));
694
695 /* xj = (*xi) + (qi - qj)/ (2.0*mi) */
696 SCIPquadprecSumDD(xjq, qi, -qj);
697 SCIPquadprecDivQD(xjq, xjq, 2.0 * mi);
698 SCIPquadprecSumQQ(xjq, xjq, xiq);
699 assert(EPSEQ(QUAD_TO_DBL(xiq) + (qi - qj)/ (2.0*mi), QUAD_TO_DBL(xjq), 1e-3));
700
701 /* yj = mj * (*xj) + qj */
702 SCIPquadprecProdQD(yjq, xjq, mj);
703 SCIPquadprecSumQD(yjq, yjq, qj);
704 assert(EPSEQ(mj * QUAD_TO_DBL(xjq) + qj, QUAD_TO_DBL(yjq), 1e-3));
705
706 /* ycoef = (*xi) + (qi - qj) / (4.0*mi) note that this is wrong in Locatelli 2016 */
707 SCIPquadprecSumDD(ycoefq, qi, -qj);
708 SCIPquadprecDivQD(ycoefq, ycoefq, 4.0 * mi);
709 SCIPquadprecSumQQ(ycoefq, ycoefq, xiq);
710 assert(EPSEQ(QUAD_TO_DBL(xiq) + (qi - qj) / (4.0*mi), QUAD_TO_DBL(ycoefq), 1e-3));
711
712 /* xcoef = 2.0*mi*(*xi) - mi * (*ycoef) + qi */
713 SCIPquadprecProdQD(xcoefq, xiq, 2.0 * mi);
714 SCIPquadprecProdQD(tmpq, ycoefq, -mi);
715 SCIPquadprecSumQQ(xcoefq, xcoefq, tmpq);
716 SCIPquadprecSumQD(xcoefq, xcoefq, qi);
717 assert(EPSEQ(2.0*mi*QUAD_TO_DBL(xiq) - mi * QUAD_TO_DBL(ycoefq) + qi, QUAD_TO_DBL(xcoefq), 1e-3));
718
719 /* constant = -mj*SQR(*xj) - (*ycoef) * qj */
720 SCIPquadprecSquareQ(constantq, xjq);
721 SCIPquadprecProdQD(constantq, constantq, -mj);
722 SCIPquadprecProdQD(tmpq, ycoefq, -qj);
723 SCIPquadprecSumQQ(constantq, constantq, tmpq);
724 /* assert(EPSEQ(-mj*SQR(QUAD_TO_DBL(xjq)) - QUAD_TO_DBL(ycoefq) * qj, QUAD_TO_DBL(constantq), 1e-3)); */
725
726 *xi = QUAD_TO_DBL(xiq);
727 *yi = QUAD_TO_DBL(yiq);
728 *xj = QUAD_TO_DBL(xjq);
729 *yj = QUAD_TO_DBL(yjq);
730 *ycoef = QUAD_TO_DBL(ycoefq);
731 *xcoef = QUAD_TO_DBL(xcoefq);
732 *constant = QUAD_TO_DBL(constantq);
733 }
734 else if( mi > 0.0 )
735 {
736 assert(mj > 0.0);
737
738 /* xi = (y + SQRT(mi*mj)*x - qi) / (REALABS(mi) + SQRT(mi*mj)) */
739 SCIPquadprecProdDD(xiq, mi, mj);
740 SCIPquadprecSqrtQ(xiq, xiq);
741 SCIPquadprecProdQD(xiq, xiq, x);
742 SCIPquadprecSumQD(xiq, xiq, y);
743 SCIPquadprecSumQD(xiq, xiq, -qi); /* (y + SQRT(mi*mj)*x - qi) */
744 SCIPquadprecProdDD(tmpq, mi, mj);
745 SCIPquadprecSqrtQ(tmpq, tmpq);
746 SCIPquadprecSumQD(tmpq, tmpq, REALABS(mi)); /* REALABS(mi) + SQRT(mi*mj) */
747 SCIPquadprecDivQQ(xiq, xiq, tmpq);
748 assert(EPSEQ((y + SQRT(mi*mj)*x - qi) / (REALABS(mi) + SQRT(mi*mj)), QUAD_TO_DBL(xiq), 1e-3));
749
750 /* yi = mi*(*xi) + qi */
751 SCIPquadprecProdQD(yiq, xiq, mi);
752 SCIPquadprecSumQD(yiq, yiq, qi);
753 assert(EPSEQ(mi*(QUAD_TO_DBL(xiq)) + qi, QUAD_TO_DBL(yiq), 1e-3));
754
755 /* xj = (y + SQRT(mi*mj)*x - qj) / (REALABS(mj) + SQRT(mi*mj)) */
756 SCIPquadprecProdDD(xjq, mi, mj);
757 SCIPquadprecSqrtQ(xjq, xjq);
758 SCIPquadprecProdQD(xjq, xjq, x);
759 SCIPquadprecSumQD(xjq, xjq, y);
760 SCIPquadprecSumQD(xjq, xjq, -qj); /* (y + SQRT(mi*mj)*x - qj) */
761 SCIPquadprecProdDD(tmpq, mi, mj);
762 SCIPquadprecSqrtQ(tmpq, tmpq);
763 SCIPquadprecSumQD(tmpq, tmpq, REALABS(mj)); /* REALABS(mj) + SQRT(mi*mj) */
764 SCIPquadprecDivQQ(xjq, xjq, tmpq);
765 assert(EPSEQ((y + SQRT(mi*mj)*x - qj) / (REALABS(mj) + SQRT(mi*mj)), QUAD_TO_DBL(xjq), 1e-3));
766
767 /* yj = mj*(*xj) + qj */
768 SCIPquadprecProdQD(yjq, xjq, mj);
769 SCIPquadprecSumQD(yjq, yjq, qj);
770 assert(EPSEQ(mj*QUAD_TO_DBL(xjq) + qj, QUAD_TO_DBL(yjq), 1e-3));
771
772 /* ycoef = (2.0*mj*(*xj) + qj - 2.0*mi*(*xi) - qi) / (mj - mi) */
773 SCIPquadprecProdQD(ycoefq, xjq, 2.0 * mj);
774 SCIPquadprecSumQD(ycoefq, ycoefq, qj);
775 SCIPquadprecProdQD(tmpq, xiq, -2.0 * mi);
776 SCIPquadprecSumQQ(ycoefq, ycoefq, tmpq);
777 SCIPquadprecSumQD(ycoefq, ycoefq, -qi);
778 SCIPquadprecSumDD(tmpq, mj, -mi);
779 SCIPquadprecDivQQ(ycoefq, ycoefq, tmpq);
780 assert(EPSEQ((2.0*mj*QUAD_TO_DBL(xjq) + qj - 2.0*mi*QUAD_TO_DBL(xiq) - qi) / (mj - mi), QUAD_TO_DBL(ycoefq), 1e-3));
781
782 /* xcoef = 2.0*mj*(*xj) + qj - mj*(*ycoef) */
783 SCIPquadprecProdQD(xcoefq, xjq, 2.0 * mj);
784 SCIPquadprecSumQD(xcoefq, xcoefq, qj);
785 SCIPquadprecProdQD(tmpq, ycoefq, -mj);
786 SCIPquadprecSumQQ(xcoefq, xcoefq, tmpq);
787 assert(EPSEQ(2.0*mj*QUAD_TO_DBL(xjq) + qj - mj*QUAD_TO_DBL(ycoefq), QUAD_TO_DBL(xcoefq), 1e-3));
788
789 /* constant = -mj*SQR(*xj) - (*ycoef) * qj */
790 SCIPquadprecSquareQ(constantq, xjq);
791 SCIPquadprecProdQD(constantq, constantq, -mj);
792 SCIPquadprecProdQD(tmpq, ycoefq, -qj);
793 SCIPquadprecSumQQ(constantq, constantq, tmpq);
794 /* assert(EPSEQ(-mj*SQR(QUAD_TO_DBL(xjq)) - QUAD_TO_DBL(ycoefq) * qj, QUAD_TO_DBL(constantq), 1e-3)); */
795
796 *xi = QUAD_TO_DBL(xiq);
797 *yi = QUAD_TO_DBL(yiq);
798 *xj = QUAD_TO_DBL(xjq);
799 *yj = QUAD_TO_DBL(yjq);
800 *ycoef = QUAD_TO_DBL(ycoefq);
801 *xcoef = QUAD_TO_DBL(xcoefq);
802 *constant = QUAD_TO_DBL(constantq);
803 }
804 else
805 {
806 assert(mi < 0.0 && mj < 0.0);
807
808 /* apply variable transformation x = -x in case for overestimation */
809 computeBilinEnvelope2(scip, -x, y, -mi, qi, -mj, qj, xi, yi, xj, yj, xcoef, ycoef, constant);
810
811 /* revert transformation; multiply cut by -1 and change -x by x */
812 *xi = -(*xi);
813 *xj = -(*xj);
814 *ycoef = -(*ycoef);
815 *constant = -(*constant);
816 }
817 }
818
819 /** computes coefficients of linearization of a bilinear term in a reference point when given two linear inequality
820 * involving only the variables of the bilinear term
821 *
822 * @note the formulas are extracted from "Convex envelopes of bivariate functions through the solution of KKT systems"
823 * by Marco Locatelli
824 *
825 */
SCIPcomputeBilinEnvelope2(SCIP * scip,SCIP_Real bilincoef,SCIP_Real lbx,SCIP_Real ubx,SCIP_Real refpointx,SCIP_Real lby,SCIP_Real uby,SCIP_Real refpointy,SCIP_Bool overestimate,SCIP_Real xcoef1,SCIP_Real ycoef1,SCIP_Real constant1,SCIP_Real xcoef2,SCIP_Real ycoef2,SCIP_Real constant2,SCIP_Real * RESTRICT lincoefx,SCIP_Real * RESTRICT lincoefy,SCIP_Real * RESTRICT linconstant,SCIP_Bool * RESTRICT success)826 void SCIPcomputeBilinEnvelope2(
827 SCIP* scip, /**< SCIP data structure */
828 SCIP_Real bilincoef, /**< coefficient of bilinear term */
829 SCIP_Real lbx, /**< lower bound on first variable */
830 SCIP_Real ubx, /**< upper bound on first variable */
831 SCIP_Real refpointx, /**< reference point for first variable */
832 SCIP_Real lby, /**< lower bound on second variable */
833 SCIP_Real uby, /**< upper bound on second variable */
834 SCIP_Real refpointy, /**< reference point for second variable */
835 SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
836 SCIP_Real xcoef1, /**< x coefficient of linear inequality; must be in {-1,0,1} */
837 SCIP_Real ycoef1, /**< y coefficient of linear inequality */
838 SCIP_Real constant1, /**< constant of linear inequality */
839 SCIP_Real xcoef2, /**< x coefficient of linear inequality; must be in {-1,0,1} */
840 SCIP_Real ycoef2, /**< y coefficient of linear inequality */
841 SCIP_Real constant2, /**< constant of linear inequality */
842 SCIP_Real* RESTRICT lincoefx, /**< buffer to store coefficient of first variable in linearization */
843 SCIP_Real* RESTRICT lincoefy, /**< buffer to store coefficient of second variable in linearization */
844 SCIP_Real* RESTRICT linconstant, /**< buffer to store constant of linearization */
845 SCIP_Bool* RESTRICT success /**< buffer to store whether linearization was successful */
846 )
847 {
848 SCIP_Real mi, mj, qi, qj, xi, xj, yi, yj;
849 SCIP_Real xcoef, ycoef, constant;
850 SCIP_Real minx, maxx, miny, maxy;
851
852 assert(scip != NULL);
853 assert(!SCIPisInfinity(scip, lbx));
854 assert(!SCIPisInfinity(scip, -ubx));
855 assert(!SCIPisInfinity(scip, lby));
856 assert(!SCIPisInfinity(scip, -uby));
857 assert(SCIPisLE(scip, lbx, ubx));
858 assert(SCIPisLE(scip, lby, uby));
859 assert(SCIPisLE(scip, lbx, refpointx));
860 assert(SCIPisGE(scip, ubx, refpointx));
861 assert(SCIPisLE(scip, lby, refpointy));
862 assert(SCIPisGE(scip, uby, refpointy));
863 assert(lincoefx != NULL);
864 assert(lincoefy != NULL);
865 assert(linconstant != NULL);
866 assert(success != NULL);
867 assert(xcoef1 != 0.0 && xcoef1 != SCIP_INVALID); /*lint !e777*/
868 assert(ycoef1 != SCIP_INVALID && ycoef1 != 0.0); /*lint !e777*/
869 assert(constant1 != SCIP_INVALID); /*lint !e777*/
870 assert(xcoef2 != 0.0 && xcoef2 != SCIP_INVALID); /*lint !e777*/
871 assert(ycoef2 != SCIP_INVALID && ycoef2 != 0.0); /*lint !e777*/
872 assert(constant2 != SCIP_INVALID); /*lint !e777*/
873
874 *success = FALSE;
875 *lincoefx = SCIP_INVALID;
876 *lincoefy = SCIP_INVALID;
877 *linconstant = SCIP_INVALID;
878
879 /* reference point does not satisfy linear inequalities */
880 if( SCIPisFeasGT(scip, xcoef1 * refpointx - ycoef1 * refpointy - constant1, 0.0)
881 || SCIPisFeasGT(scip, xcoef2 * refpointx - ycoef2 * refpointy - constant2, 0.0) )
882 return;
883
884 /* compute minimal and maximal bounds on x and y for accepting the reference point */
885 minx = lbx + 0.01 * (ubx-lbx);
886 maxx = ubx - 0.01 * (ubx-lbx);
887 miny = lby + 0.01 * (uby-lby);
888 maxy = uby - 0.01 * (uby-lby);
889
890 /* check the reference point is in the interior of the domain */
891 if( SCIPisLE(scip, refpointx, minx) || SCIPisGE(scip, refpointx, maxx)
892 || SCIPisLE(scip, refpointy, miny) || SCIPisFeasGE(scip, refpointy, maxy) )
893 return;
894
895 /* the sign of the x-coefficients of the two inequalities must be different; otherwise the convex or concave
896 * envelope can be computed via SCIPcomputeBilinEnvelope1 for each inequality separately
897 */
898 if( (xcoef1 > 0) == (xcoef2 > 0) )
899 return;
900
901 /* always consider xy without the bilinear coefficient */
902 if( bilincoef < 0.0 )
903 overestimate = !overestimate;
904
905 /* we use same notation as in "Convex envelopes of bivariate functions through the solution of KKT systems", 2016 */
906 mi = xcoef1 / ycoef1;
907 qi = -constant1 / ycoef1;
908 mj = xcoef2 / ycoef2;
909 qj = -constant2 / ycoef2;
910
911 /* mi, mj > 0 => underestimate; mi, mj < 0 => overestimate */
912 if( SCIPisNegative(scip, mi) != overestimate || SCIPisNegative(scip, mj) != overestimate )
913 return;
914
915 /* compute cut according to Locatelli 2016 */
916 computeBilinEnvelope2(scip, refpointx, refpointy, mi, qi, mj, qj, &xi, &yi, &xj, &yj, &xcoef, &ycoef, &constant);
917 assert(SCIPisRelEQ(scip, mi*xi + qi, yi));
918 assert(SCIPisRelEQ(scip, mj*xj + qj, yj));
919
920 /* it might happen that (xi,yi) = (xj,yj) if the two lines intersect */
921 if( SCIPisEQ(scip, xi, xj) && SCIPisEQ(scip, yi, yj) )
922 return;
923
924 /* check whether projected points are in the interior */
925 if( SCIPisLE(scip, xi, minx) || SCIPisGE(scip, xi, maxx) || SCIPisLE(scip, yi, miny) || SCIPisGE(scip, yi, maxy) )
926 return;
927 if( SCIPisLE(scip, xj, minx) || SCIPisGE(scip, xj, maxx) || SCIPisLE(scip, yj, miny) || SCIPisGE(scip, yj, maxy) )
928 return;
929
930 *lincoefx = bilincoef * xcoef;
931 *lincoefy = bilincoef * ycoef;
932 *linconstant = bilincoef * constant;
933
934 /* cut needs to be tight at (vx,vy) and (xj,yj) */
935 *success = SCIPisFeasEQ(scip, (*lincoefx)*xi + (*lincoefy)*yi + (*linconstant), bilincoef*xi*yi)
936 && SCIPisFeasEQ(scip, (*lincoefx)*xj + (*lincoefy)*yj + (*linconstant), bilincoef*xj*yj);
937
938 #ifndef NDEBUG
939 {
940 SCIP_Real activity = (*lincoefx)*refpointx + (*lincoefy)*refpointy + (*linconstant);
941
942 /* cut needs to under- or overestimate the bilinear term at the reference point */
943 if( bilincoef < 0.0 )
944 overestimate = !overestimate;
945
946 if( overestimate )
947 assert(SCIPisFeasGE(scip, activity, bilincoef*refpointx*refpointy));
948 else
949 assert(SCIPisFeasLE(scip, activity, bilincoef*refpointx*refpointy));
950 }
951 #endif
952 }
953
954 /** creates an NLP relaxation and stores it in a given NLPI problem; the function computes for each variable which the
955 * number of non-linearly occurrences and stores it in the nlscore array
956 *
957 * @note the first row corresponds always to the cutoff row (even if cutoffbound is SCIPinfinity(scip))
958 **/
SCIPcreateNlpiProb(SCIP * scip,SCIP_NLPI * nlpi,SCIP_NLROW ** nlrows,int nnlrows,SCIP_NLPIPROBLEM * nlpiprob,SCIP_HASHMAP * var2idx,SCIP_HASHMAP * nlrow2idx,SCIP_Real * nlscore,SCIP_Real cutoffbound,SCIP_Bool setobj,SCIP_Bool onlyconvex)959 SCIP_RETCODE SCIPcreateNlpiProb(
960 SCIP* scip, /**< SCIP data structure */
961 SCIP_NLPI* nlpi, /**< interface to NLP solver */
962 SCIP_NLROW** nlrows, /**< nonlinear rows */
963 int nnlrows, /**< total number of nonlinear rows */
964 SCIP_NLPIPROBLEM* nlpiprob, /**< empty nlpi problem */
965 SCIP_HASHMAP* var2idx, /**< empty hash map to store mapping between variables and indices in nlpi
966 * problem */
967 SCIP_HASHMAP* nlrow2idx, /**< empty hash map to store mapping between variables and indices in nlpi
968 * problem, can be NULL */
969 SCIP_Real* nlscore, /**< array to store the score of each nonlinear variable (NULL if not
970 * needed) */
971 SCIP_Real cutoffbound, /**< cutoff bound */
972 SCIP_Bool setobj, /**< should the objective function be set? */
973 SCIP_Bool onlyconvex /**< filter only for convex constraints */
974 )
975 {
976 SCIP_EXPRTREE** exprtrees;
977 int** exprvaridxs;
978 SCIP_QUADELEM** quadelems;
979 int* nquadelems;
980 SCIP_Real** linvals;
981 int** lininds;
982 int* nlininds;
983 SCIP_Real* lhss;
984 SCIP_Real* rhss;
985 const char** names;
986 SCIP_VAR** vars;
987 int nvars;
988 SCIP_Real* lbs;
989 SCIP_Real* ubs;
990 SCIP_Real* objvals = NULL;
991 int* objinds = NULL;
992 const char** varnames;
993 int nobjinds;
994 int nconss;
995 int i;
996
997 assert(nlpiprob != NULL);
998 assert(var2idx != NULL);
999 assert(nlrows != NULL);
1000 assert(nnlrows > 0);
1001 assert(nlpi != NULL);
1002
1003 SCIPdebugMsg(scip, "call SCIPcreateConvexNlpNlobbt() with cutoffbound %g\n", cutoffbound);
1004
1005 if( nlscore != NULL )
1006 {
1007 BMSclearMemoryArray(nlscore, SCIPgetNVars(scip));
1008 }
1009 vars = SCIPgetVars(scip);
1010 nvars = SCIPgetNVars(scip);
1011 nconss = 0;
1012
1013 SCIP_CALL( SCIPallocBufferArray(scip, &exprtrees, nnlrows + 1) );
1014 SCIP_CALL( SCIPallocBufferArray(scip, &exprvaridxs, nnlrows + 1) );
1015 SCIP_CALL( SCIPallocBufferArray(scip, &quadelems, nnlrows + 1) );
1016 SCIP_CALL( SCIPallocBufferArray(scip, &nquadelems, nnlrows + 1) );
1017 SCIP_CALL( SCIPallocBufferArray(scip, &linvals, nnlrows + 1) );
1018 SCIP_CALL( SCIPallocBufferArray(scip, &lininds, nnlrows + 1) );
1019 SCIP_CALL( SCIPallocBufferArray(scip, &nlininds, nnlrows + 1) );
1020 SCIP_CALL( SCIPallocBufferArray(scip, &names, nnlrows + 1) );
1021 SCIP_CALL( SCIPallocBufferArray(scip, &lhss, nnlrows + 1) );
1022 SCIP_CALL( SCIPallocBufferArray(scip, &rhss, nnlrows + 1) );
1023
1024 if( setobj )
1025 {
1026 SCIP_CALL( SCIPallocBufferArray(scip, &objvals, nvars) );
1027 SCIP_CALL( SCIPallocBufferArray(scip, &objinds, nvars) );
1028 }
1029
1030 SCIP_CALL( SCIPallocBufferArray(scip, &lbs, nvars) );
1031 SCIP_CALL( SCIPallocBufferArray(scip, &ubs, nvars) );
1032 SCIP_CALL( SCIPallocBufferArray(scip, &varnames, nvars) );
1033
1034 /* create a unique mapping between variables and {0,..,nvars-1} */
1035 nobjinds = 0;
1036 for( i = 0; i < nvars; ++i )
1037 {
1038 assert(vars[i] != NULL);
1039 SCIP_CALL( SCIPhashmapInsertInt(var2idx, (void*)vars[i], i) );
1040
1041 lbs[i] = SCIPvarGetLbLocal(vars[i]);
1042 ubs[i] = SCIPvarGetUbLocal(vars[i]);
1043 varnames[i] = SCIPvarGetName(vars[i]);
1044
1045 /* collect non-zero objective coefficients */
1046 if( setobj && !SCIPisZero(scip, SCIPvarGetObj(vars[i])) )
1047 {
1048 assert(objvals != NULL);
1049 assert(objinds != NULL);
1050
1051 objvals[nobjinds] = SCIPvarGetObj(vars[i]);
1052 objinds[nobjinds] = i;
1053 ++nobjinds;
1054 }
1055 }
1056
1057 /* add variables */
1058 SCIP_CALL( SCIPnlpiAddVars(nlpi, nlpiprob, nvars, lbs, ubs, varnames) );
1059 SCIPfreeBufferArray(scip, &varnames);
1060 SCIPfreeBufferArray(scip, &ubs);
1061 SCIPfreeBufferArray(scip, &lbs);
1062
1063 /* set the objective function */
1064 if( setobj )
1065 {
1066 if( nobjinds > 0 )
1067 {
1068 SCIP_CALL( SCIPnlpiSetObjective(nlpi, nlpiprob, nobjinds, objinds, objvals, 0, NULL, NULL, NULL, 0.0) );
1069 }
1070
1071 SCIPfreeBufferArray(scip, &objinds);
1072 SCIPfreeBufferArray(scip, &objvals);
1073 }
1074
1075 /* add row for cutoff bound even if cutoffbound == SCIPinfinity() */
1076 lhss[nconss] = -SCIPinfinity(scip);
1077 rhss[nconss] = cutoffbound;
1078 names[nconss] = "objcutoff";
1079 lininds[nconss] = NULL;
1080 linvals[nconss] = NULL;
1081 nlininds[nconss] = 0;
1082 nquadelems[nconss] = 0;
1083 quadelems[nconss] = NULL;
1084 exprtrees[nconss] = NULL;
1085 exprvaridxs[nconss] = NULL;
1086
1087 SCIP_CALL( SCIPallocBufferArray(scip, &lininds[nconss], nvars) ); /*lint !e866*/
1088 SCIP_CALL( SCIPallocBufferArray(scip, &linvals[nconss], nvars) ); /*lint !e866*/
1089
1090 for( i = 0; i < nvars; ++i )
1091 {
1092 if( !SCIPisZero(scip, SCIPvarGetObj(vars[i])) )
1093 {
1094 linvals[nconss][nlininds[nconss]] = SCIPvarGetObj(vars[i]);
1095 lininds[nconss][nlininds[nconss]] = i;
1096 ++nlininds[nconss];
1097 }
1098 }
1099 ++nconss;
1100
1101 /* add convex nonlinear rows to NLPI problem */
1102 for( i = 0; i < nnlrows; ++i )
1103 {
1104 SCIP_Bool userhs;
1105 SCIP_Bool uselhs;
1106 int k;
1107 SCIP_NLROW* nlrow;
1108
1109 nlrow = nlrows[i];
1110 assert(nlrow != NULL);
1111
1112 uselhs = FALSE;
1113 userhs = FALSE;
1114
1115 /* check curvature together with constraint sides of a nonlinear row */
1116 if( SCIPnlrowGetNQuadElems(nlrow) == 0 && SCIPnlrowGetExprtree(nlrow) == NULL )
1117 {
1118 uselhs = TRUE;
1119 userhs = TRUE;
1120 }
1121 else
1122 {
1123 if( (!onlyconvex || SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONVEX)
1124 && !SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow)) )
1125 userhs = TRUE;
1126 if( (!onlyconvex || SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONCAVE)
1127 && !SCIPisInfinity(scip, SCIPnlrowGetLhs(nlrow)) )
1128 uselhs = TRUE;
1129 }
1130
1131 if( !uselhs && !userhs )
1132 continue;
1133
1134 lhss[nconss] = uselhs ? SCIPnlrowGetLhs(nlrow) - SCIPnlrowGetConstant(nlrow) : -SCIPinfinity(scip);
1135 rhss[nconss] = userhs ? SCIPnlrowGetRhs(nlrow) - SCIPnlrowGetConstant(nlrow) : SCIPinfinity(scip);
1136 names[nconss] = SCIPnlrowGetName(nlrow);
1137 nlininds[nconss] = 0;
1138 lininds[nconss] = NULL;
1139 linvals[nconss] = NULL;
1140 nquadelems[nconss] = 0;
1141 quadelems[nconss] = NULL;
1142 exprtrees[nconss] = NULL;
1143 exprvaridxs[nconss] = NULL;
1144
1145 /* copy linear part */
1146 if( SCIPnlrowGetNLinearVars(nlrow) > 0 )
1147 {
1148 SCIP_VAR* var;
1149
1150 nlininds[nconss] = SCIPnlrowGetNLinearVars(nlrow);
1151
1152 SCIP_CALL( SCIPallocBufferArray(scip, &lininds[nconss], nlininds[nconss]) ); /*lint !e866*/
1153 SCIP_CALL( SCIPallocBufferArray(scip, &linvals[nconss], nlininds[nconss]) ); /*lint !e866*/
1154
1155 for( k = 0; k < nlininds[nconss]; ++k )
1156 {
1157 var = SCIPnlrowGetLinearVars(nlrow)[k];
1158 assert(var != NULL);
1159 assert(SCIPhashmapExists(var2idx, (void*)var));
1160
1161 lininds[nconss][k] = SCIPhashmapGetImageInt(var2idx, (void*)var);
1162 assert(var == vars[lininds[nconss][k]]);
1163 linvals[nconss][k] = SCIPnlrowGetLinearCoefs(nlrow)[k];
1164 }
1165 }
1166
1167 /* copy quadratic part */
1168 if( SCIPnlrowGetNQuadElems(nlrow) > 0 )
1169 {
1170 SCIP_QUADELEM quadelem;
1171 SCIP_VAR* var1;
1172 SCIP_VAR* var2;
1173
1174 nquadelems[nconss] = SCIPnlrowGetNQuadElems(nlrow);
1175 SCIP_CALL( SCIPallocBufferArray(scip, &quadelems[nconss], nquadelems[nconss]) ); /*lint !e866*/
1176
1177 for( k = 0; k < nquadelems[nconss]; ++k )
1178 {
1179 quadelem = SCIPnlrowGetQuadElems(nlrow)[k];
1180
1181 var1 = SCIPnlrowGetQuadVars(nlrow)[quadelem.idx1];
1182 assert(var1 != NULL);
1183 assert(SCIPhashmapExists(var2idx, (void*)var1));
1184
1185 var2 = SCIPnlrowGetQuadVars(nlrow)[quadelem.idx2];
1186 assert(var2 != NULL);
1187 assert(SCIPhashmapExists(var2idx, (void*)var2));
1188
1189 quadelems[nconss][k].coef = quadelem.coef;
1190 quadelems[nconss][k].idx1 = SCIPhashmapGetImageInt(var2idx, (void*)var1);
1191 quadelems[nconss][k].idx2 = SCIPhashmapGetImageInt(var2idx, (void*)var2);
1192
1193 /* expr.c assumes that the indices are ordered */
1194 if( quadelems[nconss][k].idx1 > quadelems[nconss][k].idx2 )
1195 {
1196 SCIPswapInts(&quadelems[nconss][k].idx1, &quadelems[nconss][k].idx2);
1197 }
1198 assert(quadelems[nconss][k].idx1 <= quadelems[nconss][k].idx2);
1199
1200 /* update nlscore */
1201 if( nlscore != NULL )
1202 {
1203 ++nlscore[quadelems[nconss][k].idx1];
1204 if( quadelems[nconss][k].idx1 != quadelems[nconss][k].idx2 )
1205 ++nlscore[quadelems[nconss][k].idx2];
1206 }
1207 }
1208 }
1209
1210 /* copy expression tree */
1211 if( SCIPnlrowGetExprtree(nlrow) != NULL )
1212 {
1213 SCIP_VAR* var;
1214
1215 /* note that we don't need to copy the expression tree here since only the mapping between variables in the
1216 * tree and the corresponding indices change; this mapping is stored in the exprvaridxs array
1217 */
1218 exprtrees[nconss] = SCIPnlrowGetExprtree(nlrow);
1219
1220 SCIP_CALL( SCIPallocBufferArray(scip, &exprvaridxs[nconss], SCIPexprtreeGetNVars(exprtrees[nconss])) ); /*lint !e866*/
1221
1222 for( k = 0; k < SCIPexprtreeGetNVars(exprtrees[nconss]); ++k )
1223 {
1224 var = SCIPexprtreeGetVars(exprtrees[nconss])[k];
1225 assert(var != NULL);
1226 assert(SCIPhashmapExists(var2idx, (void*)var));
1227
1228 exprvaridxs[nconss][k] = SCIPhashmapGetImageInt(var2idx, (void*)var);
1229
1230 /* update nlscore */
1231 if( nlscore != NULL )
1232 ++nlscore[exprvaridxs[nconss][k]];
1233 }
1234 }
1235
1236 /* if the row to index hash map is provided, we need to store the row index */
1237 if( nlrow2idx != NULL )
1238 {
1239 SCIP_CALL( SCIPhashmapInsertInt(nlrow2idx, nlrow, nconss) );
1240 }
1241
1242 ++nconss;
1243 }
1244 assert(nconss > 0);
1245
1246 /* pass all constraint information to nlpi */
1247 SCIP_CALL( SCIPnlpiAddConstraints(nlpi, nlpiprob, nconss, lhss, rhss, nlininds, lininds, linvals, nquadelems,
1248 quadelems, exprvaridxs, exprtrees, names) );
1249
1250 /* free memory */
1251 for( i = nconss - 1; i > 0; --i )
1252 {
1253 if( exprtrees[i] != NULL )
1254 {
1255 assert(exprvaridxs[i] != NULL);
1256 SCIPfreeBufferArray(scip, &exprvaridxs[i]);
1257 }
1258
1259 if( nquadelems[i] > 0 )
1260 {
1261 assert(quadelems[i] != NULL);
1262 SCIPfreeBufferArray(scip, &quadelems[i]);
1263 }
1264
1265 if( nlininds[i] > 0 )
1266 {
1267 assert(linvals[i] != NULL);
1268 assert(lininds[i] != NULL);
1269 SCIPfreeBufferArray(scip, &linvals[i]);
1270 SCIPfreeBufferArray(scip, &lininds[i]);
1271 }
1272 }
1273 /* free row for cutoff bound even if objective is 0 */
1274 SCIPfreeBufferArray(scip, &linvals[i]);
1275 SCIPfreeBufferArray(scip, &lininds[i]);
1276
1277 SCIPfreeBufferArray(scip, &rhss);
1278 SCIPfreeBufferArray(scip, &lhss);
1279 SCIPfreeBufferArray(scip, &names);
1280 SCIPfreeBufferArray(scip, &nlininds);
1281 SCIPfreeBufferArray(scip, &lininds);
1282 SCIPfreeBufferArray(scip, &linvals);
1283 SCIPfreeBufferArray(scip, &nquadelems);
1284 SCIPfreeBufferArray(scip, &quadelems);
1285 SCIPfreeBufferArray(scip, &exprvaridxs);
1286 SCIPfreeBufferArray(scip, &exprtrees);
1287
1288 return SCIP_OKAY;
1289 }
1290
1291 /** updates bounds of each variable and the cutoff row in the nlpiproblem */
SCIPupdateNlpiProb(SCIP * scip,SCIP_NLPI * nlpi,SCIP_NLPIPROBLEM * nlpiprob,SCIP_HASHMAP * var2nlpiidx,SCIP_VAR ** nlpivars,int nlpinvars,SCIP_Real cutoffbound)1292 SCIP_RETCODE SCIPupdateNlpiProb(
1293 SCIP* scip, /**< SCIP data structure */
1294 SCIP_NLPI* nlpi, /**< interface to NLP solver */
1295 SCIP_NLPIPROBLEM* nlpiprob, /**< nlpi problem representing the convex NLP relaxation */
1296 SCIP_HASHMAP* var2nlpiidx, /**< mapping between variables and nlpi indices */
1297 SCIP_VAR** nlpivars, /**< array containing all variables of the nlpi */
1298 int nlpinvars, /**< total number of nlpi variables */
1299 SCIP_Real cutoffbound /**< new cutoff bound */
1300 )
1301 {
1302 SCIP_Real* lbs;
1303 SCIP_Real* ubs;
1304 SCIP_Real lhs;
1305 SCIP_Real rhs;
1306 int* inds;
1307 int i;
1308
1309 SCIPdebugMsg(scip, "call SCIPupdateConvexNlpNlobbt()\n");
1310
1311 /* update variable bounds */
1312 SCIP_CALL( SCIPallocBufferArray(scip, &lbs, nlpinvars) );
1313 SCIP_CALL( SCIPallocBufferArray(scip, &ubs, nlpinvars) );
1314 SCIP_CALL( SCIPallocBufferArray(scip, &inds, nlpinvars) );
1315
1316 for( i = 0; i < nlpinvars; ++i )
1317 {
1318 assert(nlpivars[i] != NULL);
1319 assert(SCIPhashmapExists(var2nlpiidx, (void*)nlpivars[i]));
1320
1321 lbs[i] = SCIPvarGetLbLocal(nlpivars[i]);
1322 ubs[i] = SCIPvarGetUbLocal(nlpivars[i]);
1323 inds[i] = SCIPhashmapGetImageInt(var2nlpiidx, (void*)nlpivars[i]);
1324 assert(inds[i] >= 0 && inds[i] < nlpinvars);
1325 }
1326
1327 SCIP_CALL( SCIPnlpiChgVarBounds(nlpi, nlpiprob, nlpinvars, inds, lbs, ubs) );
1328
1329 SCIPfreeBufferArray(scip, &inds);
1330 SCIPfreeBufferArray(scip, &ubs);
1331 SCIPfreeBufferArray(scip, &lbs);
1332
1333 /* update cutoff row */
1334 lhs = -SCIPinfinity(scip);
1335 rhs = cutoffbound;
1336 i = 0;
1337
1338 SCIP_CALL( SCIPnlpiChgConsSides(nlpi, nlpiprob, 1, &i, &lhs, &rhs) );
1339
1340 return SCIP_OKAY;
1341 }
1342
1343 /** adds linear rows to the NLP relaxation */
SCIPaddNlpiProbRows(SCIP * scip,SCIP_NLPI * nlpi,SCIP_NLPIPROBLEM * nlpiprob,SCIP_HASHMAP * var2idx,SCIP_ROW ** rows,int nrows)1344 SCIP_RETCODE SCIPaddNlpiProbRows(
1345 SCIP* scip, /**< SCIP data structure */
1346 SCIP_NLPI* nlpi, /**< interface to NLP solver */
1347 SCIP_NLPIPROBLEM* nlpiprob, /**< nlpi problem */
1348 SCIP_HASHMAP* var2idx, /**< empty hash map to store mapping between variables and indices in nlpi
1349 * problem */
1350 SCIP_ROW** rows, /**< rows to add */
1351 int nrows /**< total number of rows to add */
1352 )
1353 {
1354 const char** names;
1355 SCIP_Real* lhss;
1356 SCIP_Real* rhss;
1357 SCIP_Real** linvals;
1358 int** lininds;
1359 int* nlininds;
1360 int i;
1361
1362 assert(nlpi != NULL);
1363 assert(nlpiprob != NULL);
1364 assert(var2idx != NULL);
1365 assert(nrows == 0 || rows != NULL);
1366
1367 SCIPdebugMsg(scip, "call SCIPaddConvexNlpRowsNlobbt() with %d rows\n", nrows);
1368
1369 if( nrows <= 0 )
1370 return SCIP_OKAY;
1371
1372 SCIP_CALL( SCIPallocBufferArray(scip, &names, nrows) );
1373 SCIP_CALL( SCIPallocBufferArray(scip, &lhss, nrows) );
1374 SCIP_CALL( SCIPallocBufferArray(scip, &rhss, nrows) );
1375 SCIP_CALL( SCIPallocBufferArray(scip, &linvals, nrows) );
1376 SCIP_CALL( SCIPallocBufferArray(scip, &lininds, nrows) );
1377 SCIP_CALL( SCIPallocBufferArray(scip, &nlininds, nrows) );
1378
1379 for( i = 0; i < nrows; ++i )
1380 {
1381 int k;
1382
1383 assert(rows[i] != NULL);
1384 assert(SCIProwGetNNonz(rows[i]) <= SCIPgetNVars(scip));
1385
1386 names[i] = SCIProwGetName(rows[i]);
1387 lhss[i] = SCIProwGetLhs(rows[i]) - SCIProwGetConstant(rows[i]);
1388 rhss[i] = SCIProwGetRhs(rows[i]) - SCIProwGetConstant(rows[i]);
1389 nlininds[i] = SCIProwGetNNonz(rows[i]);
1390 linvals[i] = SCIProwGetVals(rows[i]);
1391 lininds[i] = NULL;
1392
1393 SCIP_CALL( SCIPallocBufferArray(scip, &lininds[i], SCIProwGetNNonz(rows[i])) ); /*lint !e866*/
1394
1395 for( k = 0; k < SCIProwGetNNonz(rows[i]); ++k )
1396 {
1397 SCIP_VAR* var;
1398
1399 var = SCIPcolGetVar(SCIProwGetCols(rows[i])[k]);
1400 assert(var != NULL);
1401 assert(SCIPhashmapExists(var2idx, (void*)var));
1402
1403 lininds[i][k] = SCIPhashmapGetImageInt(var2idx, (void*)var);
1404 assert(lininds[i][k] >= 0 && lininds[i][k] < SCIPgetNVars(scip));
1405 }
1406 }
1407
1408 /* pass all linear rows to the nlpi */
1409 SCIP_CALL( SCIPnlpiAddConstraints(nlpi, nlpiprob, nrows, lhss, rhss, nlininds, lininds, linvals, NULL,
1410 NULL, NULL, NULL, names) );
1411
1412 /* free memory */
1413 for( i = nrows - 1; i >= 0; --i )
1414 {
1415 SCIPfreeBufferArray(scip, &lininds[i]);
1416 }
1417 SCIPfreeBufferArray(scip, &nlininds);
1418 SCIPfreeBufferArray(scip, &lininds);
1419 SCIPfreeBufferArray(scip, &linvals);
1420 SCIPfreeBufferArray(scip, &rhss);
1421 SCIPfreeBufferArray(scip, &lhss);
1422 SCIPfreeBufferArray(scip, &names);
1423
1424 return SCIP_OKAY;
1425 }
1426