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 gauge.c
17 * @brief unit test for sepa_gauge methods
18 */
19
20 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
21
22 #include "scip/scip.h"
23 #include "nlpi/nlpi_ipopt.h"
24
25 #include "scip/sepa_gauge.c"
26
27 #include "include/scip_test.h"
28
29 #define EPS 1e-6
30
31 static SCIP* scip = NULL;
32 static SCIP_NLPI* nlpi = NULL;
33 static SCIP_NLROW* nlrow1 = NULL;
34 static SCIP_NLROW* nlrow2 = NULL;
35 static SCIP_NLROW* nlrow3 = NULL;
36 static SCIP_VAR* x;
37 static SCIP_VAR* y;
38
39 /* methods to test:
40 * computeInteriorPoint: uses SCIP's NLP to build a convex NLP relaxation, solves it and store the solution in the sepadata.
41 * findBoundaryPoint: receives the nlrows, on which sides they are convex, indices of nlrows to separate, number of indices, the interior solution, the solution to separate and the final solution
42 * generateCut: receives point where to compute gradient, exprinterpreter (?), the nlrow, the side which is convex and where to store the cut
43 */
44
45 /* create log(exp(x) + exp(y)) <= 1 */
46 static
createNlRow1(CONVEXSIDE convexside)47 void createNlRow1(CONVEXSIDE convexside)
48 {
49 SCIP_EXPRTREE* exprtree;
50 SCIP_EXPRCURV curvature;
51 SCIP_EXPR* xexpr;
52 SCIP_EXPR* yexpr;
53 SCIP_EXPR* exp_x;
54 SCIP_EXPR* exp_y;
55 SCIP_EXPR* sum_exp;
56 SCIP_EXPR* lgexp;
57 SCIP_EXPR* finalexpr;
58 SCIP_Real lhs;
59 SCIP_Real rhs;
60 SCIP_VAR* vars[2];
61
62 /* setup expression for x and y */
63 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &xexpr, SCIP_EXPR_VARIDX, 0) );
64 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &yexpr, SCIP_EXPR_VARIDX, 1) );
65
66 /* setup expression for exp(x) and exp(y) */
67 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exp_x, SCIP_EXPR_EXP, xexpr) );
68 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exp_y, SCIP_EXPR_EXP, yexpr) );
69
70 /* setup expression for exp(x) + exp(y) */
71 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &sum_exp, SCIP_EXPR_PLUS, exp_x, exp_y) );
72
73 /* setup expression for lg(sum) */
74 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &lgexp, SCIP_EXPR_LOG, sum_exp) );
75
76 /* decide curvature */
77 if( convexside == RHS )
78 {
79 curvature = SCIP_EXPRCURV_CONVEX;
80 lhs = -SCIPinfinity(scip);
81 rhs = 1.0;
82 finalexpr = lgexp;
83 }
84 else
85 {
86 SCIP_EXPR* minusone;
87 curvature = SCIP_EXPRCURV_CONCAVE;
88 lhs = -1.0;
89 rhs = SCIPinfinity(scip);
90 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &minusone, SCIP_EXPR_CONST, -1.0) );
91 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &finalexpr, SCIP_EXPR_MUL, minusone, lgexp) );
92 }
93
94 /* setup expression tree with finalexpr as root expression */
95 SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, finalexpr, 2, 0, NULL) );
96
97 /* assign SCIP variables to tree */
98 vars[0] = x;
99 vars[1] = y;
100 SCIP_CALL( SCIPexprtreeSetVars(exprtree, 2, vars) );
101
102 /* create nonlinear row for exprtree <= 1 */
103 SCIP_CALL( SCIPcreateNlRow(scip, &nlrow1, "lg(sum exp)", 0.0, 0, NULL, NULL, 0, NULL, 0, NULL, exprtree, lhs, rhs, curvature) );
104
105 /* release */
106 SCIP_CALL( SCIPexprtreeFree(&exprtree) );
107 }
108
109 /* create x^2 <= y */
110 static
createNlRow2(CONVEXSIDE convexside)111 void createNlRow2(CONVEXSIDE convexside)
112 {
113 SCIP_EXPRTREE* exprtree;
114 SCIP_EXPRCURV curvature;
115 SCIP_EXPR* xexpr;
116 SCIP_EXPR* yexpr;
117 SCIP_EXPR* sqr_x;
118 SCIP_EXPR* diff;
119 SCIP_Real lhs;
120 SCIP_Real rhs;
121 SCIP_VAR* vars[2];
122
123 /* setup expression for x and y */
124 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &xexpr, SCIP_EXPR_VARIDX, 0) );
125 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &yexpr, SCIP_EXPR_VARIDX, 1) );
126
127 /* setup expression for x^2 */
128 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &sqr_x, SCIP_EXPR_SQUARE, xexpr) );
129
130 /* decide curvature */
131 if( convexside == RHS )
132 {
133 curvature = SCIP_EXPRCURV_CONVEX;
134 lhs = -SCIPinfinity(scip);
135 rhs = 0.0;
136 /* setup expression for x^2 - y */
137 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &diff, SCIP_EXPR_MINUS, sqr_x, yexpr) );
138 }
139 else
140 {
141 curvature = SCIP_EXPRCURV_CONCAVE;
142 lhs = 0.0;
143 rhs = SCIPinfinity(scip);
144 /* setup expression for - x^2 + y */
145 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &diff, SCIP_EXPR_MINUS, yexpr, sqr_x) );
146 }
147
148 /* setup expression trees with exprs as root expression */
149 SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, diff, 2, 0, NULL) );
150
151 /* assign SCIP variables to tree */
152 vars[0] = x;
153 vars[1] = y;
154 SCIP_CALL( SCIPexprtreeSetVars(exprtree, 2, vars) );
155
156 /* create nonlinear row for exprtree <= 0 */
157 SCIP_CALL( SCIPcreateNlRow(scip, &nlrow2, "x^2 ? y", 0.0, 0, NULL, NULL, 0, NULL, 0, NULL, exprtree, lhs, rhs, curvature) );
158
159 /* release */
160 SCIP_CALL( SCIPexprtreeFree(&exprtree) );
161 }
162
163 /* 1.1*x+2.4*x^2 + 0.01*x*y + 0.3*y^2 + 0.2*log(0.5*exp(0.12*x+0.1)+2*exp(0.1*y)+0.7) <= 0.5 */
164 static
createNlRow3(CONVEXSIDE convexside)165 void createNlRow3(CONVEXSIDE convexside)
166 {
167 SCIP_EXPRTREE* exprtree;
168 SCIP_EXPRCURV curvature;
169 SCIP_EXPR* xexpr;
170 SCIP_EXPR* arg_exp_x;
171 SCIP_EXPR* arg_exp_y;
172 SCIP_EXPR* arg_log;
173 SCIP_EXPR* yexpr;
174 SCIP_EXPR* exp_x;
175 SCIP_EXPR* exp_y;
176 SCIP_EXPR* lgexp;
177 SCIP_EXPR* finalexpr;
178 SCIP_Real lhs;
179 SCIP_Real rhs;
180 SCIP_Real side;
181 SCIP_VAR* vars[2];
182
183 /* expression for x and y */
184 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &xexpr, SCIP_EXPR_VARIDX, 0) );
185 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &yexpr, SCIP_EXPR_VARIDX, 1) );
186
187 /* expression for 0.12x+0.1 and 0.1*y */
188 SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &arg_exp_x, 1, &xexpr, (SCIP_Real[]){0.12}, 0.1) );
189 SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &arg_exp_y, 1, &yexpr, (SCIP_Real[]){0.1}, 0.0) );
190
191 /* expression for exp(0.12*x+0.1) and exp(0.1*y) */
192 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exp_x, SCIP_EXPR_EXP, arg_exp_x) );
193 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exp_y, SCIP_EXPR_EXP, arg_exp_y) );
194
195 /* expression for 0.5*exp(0.12*x+0.1) + 2*exp(0.1*y) + 0.7 */
196 SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &arg_log, 2,
197 (SCIP_EXPR*[]){exp_x, exp_y}, (SCIP_Real[]){0.5, 2.0}, 0.7) );
198
199 /* expression for lg(0.5*exp(0.12*x+0.1) + 2*exp(0.1*y) + 0.7) */
200 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &lgexp, SCIP_EXPR_LOG, arg_log) );
201
202 /* decide curvature */
203 if( convexside == RHS )
204 {
205 curvature = SCIP_EXPRCURV_CONVEX;
206 side = 1.0;
207 lhs = -SCIPinfinity(scip);
208 rhs = 0.5;
209 }
210 else
211 {
212 curvature = SCIP_EXPRCURV_CONCAVE;
213 side = -1.0;
214 lhs = -0.5;
215 rhs = SCIPinfinity(scip);
216 }
217
218 SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &finalexpr, 1,
219 &lgexp, (SCIP_Real[]){side * 0.2}, 0.0) );
220
221 /* build expression tree with log as root expression */
222 SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, finalexpr, 2, 0, NULL) );
223
224 /* assign SCIP variables to tree */
225 vars[0] = x;
226 vars[1] = y;
227 SCIP_CALL( SCIPexprtreeSetVars(exprtree, 2, vars) );
228
229 SCIP_CALL( SCIPcreateNlRow(scip, &nlrow3, "complicated",
230 0.0, 1, &x, (SCIP_Real[]){side*1.1},
231 2, vars,
232 3,
233 (SCIP_QUADELEM[]){
234 {0, 0, side*2.4}, // 2.4 * x^2
235 {0, 1, side*0.01}, // 0.01*x*y
236 {1, 1, side*0.3}
237 },
238 exprtree, lhs, rhs, curvature) );
239
240
241
242 /* release */
243 SCIP_CALL( SCIPexprtreeFree(&exprtree) );
244 }
245
246 static
evaluation_setup(void)247 void evaluation_setup(void)
248 {
249 SCIP_CALL( SCIPcreate(&scip) );
250
251 /* include quadratic conshdlr (need to include nonlinear) */
252 //SCIP_CALL( SCIPincludeConshdlrNonlinear(scip) );
253 //SCIP_CALL( SCIPincludeConshdlrQuadratic(scip) );
254
255 /* include NLPI's */
256 SCIP_CALL( SCIPcreateNlpSolverIpopt(SCIPblkmem(scip), &nlpi) );
257 /* if no IPOPT available, don't run test */
258 if( nlpi == NULL )
259 return;
260
261 SCIP_CALL( SCIPincludeNlpi(scip, nlpi) );
262 SCIP_CALL( SCIPincludeExternalCodeInformation(scip, SCIPgetSolverNameIpopt(), SCIPgetSolverDescIpopt()) );
263
264 /* create a problem */
265 SCIP_CALL( SCIPcreateProbBasic(scip, "problem") );
266
267 /* change SCIP's stage to be able to create nlrows and rows */
268 SCIP_CALL( TESTscipSetStage(scip, SCIP_STAGE_SOLVING, FALSE) );
269
270 /* create and add variables */
271 SCIP_CALL( SCIPcreateVarBasic(scip, &x, "x", -2.5, 3.1, 0.0, SCIP_VARTYPE_CONTINUOUS) );
272 SCIP_CALL( SCIPcreateVarBasic(scip, &y, "y", -2.5, 3.1, 0.0, SCIP_VARTYPE_CONTINUOUS) );
273 SCIP_CALL( SCIPaddVar(scip, x) );
274 SCIP_CALL( SCIPaddVar(scip, y) );
275
276 }
277
278 static
teardown(void)279 void teardown(void)
280 {
281 /* if no IPOPT available, don't run test */
282 if( nlpi == NULL )
283 return;
284
285 SCIP_CALL( SCIPreleaseVar(scip, &x) );
286 SCIP_CALL( SCIPreleaseVar(scip, &y) );
287 if( nlrow1 != NULL )
288 {
289 SCIP_CALL( SCIPreleaseNlRow(scip, &nlrow1) );
290 }
291 if( nlrow2 != NULL )
292 {
293 SCIP_CALL( SCIPreleaseNlRow(scip, &nlrow2) );
294 }
295 if( nlrow3 != NULL )
296 {
297 SCIP_CALL( SCIPreleaseNlRow(scip, &nlrow3) );
298 }
299
300 /* free SCIP */
301 SCIP_CALL( SCIPfree(&scip) );
302
303 /* check for memory leaks */
304 cr_assert_eq(BMSgetMemoryUsed(), 0, "There is are memory leak!!");
305 }
306
307 static
evaluate_gauge(CONVEXSIDE * convexsides)308 void evaluate_gauge(CONVEXSIDE* convexsides)
309 {
310 SCIP_SOL* interior_sol;
311 SCIP_SOL* toseparate_sol;
312 SCIP_SOL* boundary_sol;
313 POSITION position;
314 SCIP_Real feas;
315 int nlrowsidx[] = {0, 1};
316 int nnlrowsidx = 2;
317 SCIP_NLROW* nlrows[2];
318
319 /* if no IPOPT available, don't run test */
320 if( nlpi == NULL )
321 return;
322
323 /* create the nl rows */
324 createNlRow1(convexsides[0]);
325 createNlRow2(convexsides[1]);
326 nlrows[0] = nlrow1;
327 nlrows[1] = nlrow2;
328
329
330 /** first point active only on nlrow1 **/
331 /* set interior point */
332 SCIP_CALL( SCIPcreateSol(scip, &interior_sol, NULL) );
333 SCIP_CALL( SCIPsetSolVal(scip, interior_sol, x, -0.5) );
334 SCIP_CALL( SCIPsetSolVal(scip, interior_sol, y, 0.5) );
335 /* set point to separate */
336 SCIP_CALL( SCIPcreateSol(scip, &toseparate_sol, NULL) );
337 SCIP_CALL( SCIPsetSolVal(scip, toseparate_sol, x, 0.75) );
338 SCIP_CALL( SCIPsetSolVal(scip, toseparate_sol, y, 0.25) );
339
340 /* find boundary point */
341 position = EXTERIOR;
342 SCIP_CALL( SCIPcreateSol(scip, &boundary_sol, NULL) );
343 SCIP_CALL( findBoundaryPoint(scip, nlrows, nlrowsidx, nnlrowsidx, convexsides, interior_sol, toseparate_sol, boundary_sol, &position) );
344 cr_expect_eq(position, BOUNDARY, "position of boundary sol is %d, expected boundary", position);
345
346 cr_expect_float_eq(SCIPgetSolVal(scip, boundary_sol, x), 0.265033442069232, EPS, "for x got %f\n", SCIPgetSolVal(scip, boundary_sol, x));
347 cr_expect_float_eq(SCIPgetSolVal(scip, boundary_sol, y), 0.346993311586154, EPS, "for y got %f\n", SCIPgetSolVal(scip, boundary_sol, y));
348
349 SCIP_CALL( SCIPgetNlRowSolFeasibility(scip, nlrow1, boundary_sol, &feas) );
350 cr_expect_float_eq(feas, 0.0, EPS, "nlrow1 is not active: feas %f\n", feas);
351 SCIP_CALL( SCIPgetNlRowSolFeasibility(scip, nlrow2, boundary_sol, &feas) );
352 cr_expect_float_neq(feas, 0.0, EPS, "nlrow2 is active: feas %f\n", feas);
353
354 /** second point active on both **/
355 /* set interior point */
356 SCIP_CALL( SCIPsetSolVal(scip, interior_sol, x, -0.8) );
357 SCIP_CALL( SCIPsetSolVal(scip, interior_sol, y, 0.75) );
358 /* set point to separate */
359 SCIP_CALL( SCIPsetSolVal(scip, toseparate_sol, x, 0.80011239) );
360 SCIP_CALL( SCIPsetSolVal(scip, toseparate_sol, y, 0.0) );
361
362 /* find boundary point */
363 position = EXTERIOR;
364 SCIP_CALL( findBoundaryPoint(scip, nlrows, nlrowsidx, nnlrowsidx, convexsides, interior_sol, toseparate_sol, boundary_sol, &position) );
365 cr_expect_eq(position, BOUNDARY, "position of boundary sol is %d, expected boundary", position);
366
367 cr_expect_float_eq(SCIPgetSolVal(scip, boundary_sol, x), 0.4213473910790077, EPS, "for x got %f\n", SCIPgetSolVal(scip, boundary_sol, x));
368 cr_expect_float_eq(SCIPgetSolVal(scip, boundary_sol, y), 0.1775336239690863, EPS, "for y got %f\n", SCIPgetSolVal(scip, boundary_sol, y));
369
370 SCIP_CALL( SCIPgetNlRowSolFeasibility(scip, nlrow1, boundary_sol, &feas) );
371 cr_expect_float_eq(feas, 0.0, EPS, "nlrow1 is not active: feas %f\n", feas);
372 SCIP_CALL( SCIPgetNlRowSolFeasibility(scip, nlrow2, boundary_sol, &feas) );
373 cr_expect_float_eq(feas, 0.0, EPS, "nlrow2 is not active: feas %f\n", feas);
374
375 /** third point active on nlrow2 **/
376 /* set interior point */
377 SCIP_CALL( SCIPsetSolVal(scip, interior_sol, x, -0.5) );
378 SCIP_CALL( SCIPsetSolVal(scip, interior_sol, y, 0.5) );
379 /* set point to separate */
380 SCIP_CALL( SCIPsetSolVal(scip, toseparate_sol, x, -1.0) );
381 SCIP_CALL( SCIPsetSolVal(scip, toseparate_sol, y, 0.0) );
382
383 /* find boundary point */
384 position = EXTERIOR;
385 SCIP_CALL( findBoundaryPoint(scip, nlrows, nlrowsidx, nnlrowsidx, convexsides, interior_sol, toseparate_sol, boundary_sol, &position) );
386 cr_expect_eq(position, BOUNDARY, "position of boundary sol is %d, expected boundary", position);
387
388 cr_expect_float_eq(SCIPgetSolVal(scip, boundary_sol, x), -0.618033988749895, EPS, "for x got %f\n", SCIPgetSolVal(scip, boundary_sol, x));
389 cr_expect_float_eq(SCIPgetSolVal(scip, boundary_sol, y), 0.381966011250105, EPS, "for y got %f\n", SCIPgetSolVal(scip, boundary_sol, y));
390
391 SCIP_CALL( SCIPgetNlRowSolFeasibility(scip, nlrow1, boundary_sol, &feas) );
392 cr_expect_float_neq(feas, 0.0, EPS, "nlrow1 is active: feas %f\n", feas);
393 SCIP_CALL( SCIPgetNlRowSolFeasibility(scip, nlrow2, boundary_sol, &feas) );
394 cr_expect_float_eq(feas, 0.0, EPS, "nlrow2 is not active: feas %f\n", feas);
395
396 SCIP_CALL( SCIPfreeSol(scip, &interior_sol) );
397 SCIP_CALL( SCIPfreeSol(scip, &toseparate_sol) );
398 SCIP_CALL( SCIPfreeSol(scip, &boundary_sol) );
399 }
400
401 /* TEST SUITE */
402 TestSuite(evaluation, .init = evaluation_setup, .fini = teardown);
403
404 /* Test: with nlrows having different convexity */
Test(evaluation,convex_convex)405 Test(evaluation, convex_convex)
406 {
407 evaluate_gauge((CONVEXSIDE[]){RHS, RHS});
408 }
Test(evaluation,concave_convex)409 Test(evaluation, concave_convex)
410 {
411 evaluate_gauge((CONVEXSIDE[]){LHS, RHS});
412 }
Test(evaluation,concave_concave)413 Test(evaluation, concave_concave)
414 {
415 evaluate_gauge((CONVEXSIDE[]){LHS, LHS});
416 }
Test(evaluation,convex_concave)417 Test(evaluation, convex_concave)
418 {
419 evaluate_gauge((CONVEXSIDE[]){RHS, LHS});
420 }
421
422 /* test that we get the correct gradient cut for convex regions, defined as both, convex <= rhs and concave >= lhs */
Test(evaluation,gradient_cut_convex)423 Test(evaluation, gradient_cut_convex)
424 {
425 SCIP_SOL* x0;
426 SCIP_ROW* gradcut;
427 SCIP_EXPRINT* exprint;
428 SCIP_Bool convex = TRUE;
429 SCIP_Real* coefs;
430 SCIP_Bool success = FALSE;
431
432 /* if no IPOPT available, don't run test */
433 if( nlpi == NULL )
434 return;
435
436 /* create the nl row */
437 createNlRow1(convex);
438
439 /** point where to compute gradient cut **/
440 SCIP_CALL( SCIPcreateSol(scip, &x0, NULL) );
441 SCIP_CALL( SCIPsetSolVal(scip, x0, x, -0.5) );
442 SCIP_CALL( SCIPsetSolVal(scip, x0, y, 0.5) );
443
444 /* compute gradient cut */
445 SCIP_CALL( SCIPcreateEmptyRowUnspec(scip, &gradcut, "gradcut", -SCIPinfinity(scip), SCIPinfinity(scip),
446 FALSE, FALSE , FALSE) );
447 SCIP_CALL( SCIPexprintCreate(SCIPblkmem(scip), &exprint) );
448 SCIP_CALL( generateCut(scip, x0, exprint, nlrow1, convex, gradcut, &success) );
449 cr_assert(success);
450
451 /* check coefficients */
452 coefs = SCIProwGetVals(gradcut);
453 cr_expect_eq(SCIProwGetNNonz(gradcut), 2);
454 cr_expect_float_eq(0.268941421369995, coefs[0], EPS, "for x got %f\n", coefs[0]);
455 cr_expect_float_eq(0.731058578630005, coefs[1], EPS, "for y got %f\n", coefs[1]);
456 cr_expect_float_eq(0.417796891111782, SCIProwGetRhs(gradcut), EPS, "wrong rhs\n");
457
458 SCIP_CALL( SCIPexprintFree(&exprint) );
459 SCIP_CALL( SCIPfreeSol(scip, &x0) );
460 SCIP_CALL( SCIPreleaseRow(scip, &gradcut) );
461 }
462
Test(evaluation,gradient_cut_concave)463 Test(evaluation, gradient_cut_concave)
464 {
465 SCIP_SOL* x0;
466 SCIP_ROW* gradcut;
467 SCIP_EXPRINT* exprint;
468 SCIP_Bool convex = FALSE;
469 SCIP_Real* coefs;
470 SCIP_Bool success = FALSE;
471
472 /* if no IPOPT available, don't run test */
473 if( nlpi == NULL )
474 return;
475
476 /* create the nl row */
477 createNlRow1(convex);
478
479 /** point where to compute gradient cut **/
480 SCIP_CALL( SCIPcreateSol(scip, &x0, NULL) );
481 SCIP_CALL( SCIPsetSolVal(scip, x0, x, -0.5) );
482 SCIP_CALL( SCIPsetSolVal(scip, x0, y, 0.5) );
483
484 /* compute gradient cut */
485 SCIP_CALL( SCIPcreateEmptyRowUnspec(scip, &gradcut, "gradcut", -SCIPinfinity(scip), SCIPinfinity(scip),
486 FALSE, FALSE , FALSE) );
487 SCIP_CALL( SCIPexprintCreate(SCIPblkmem(scip), &exprint) );
488 SCIP_CALL( generateCut(scip, x0, exprint, nlrow1, convex, gradcut, &success) );
489 cr_assert(success);
490
491 /* check coefficients */
492 coefs = SCIProwGetVals(gradcut);
493 cr_expect_eq(SCIProwGetNNonz(gradcut), 2);
494 cr_expect_float_eq(-0.268941421369995, coefs[0], EPS, "for x got %f\n", coefs[0]);
495 cr_expect_float_eq(-0.731058578630005, coefs[1], EPS, "for y got %f\n", coefs[1]);
496 cr_expect_float_eq(-0.417796891111782, SCIProwGetLhs(gradcut), EPS, "wrong rhs\n");
497
498 SCIP_CALL( SCIPexprintFree(&exprint) );
499 SCIP_CALL( SCIPfreeSol(scip, &x0) );
500 SCIP_CALL( SCIPreleaseRow(scip, &gradcut) );
501 }
502
Test(evaluation,gradient_complicated_convex)503 Test(evaluation, gradient_complicated_convex)
504 {
505 SCIP_SOL* x0;
506 SCIP_ROW* gradcut;
507 SCIP_EXPRINT* exprint;
508 SCIP_Bool convex = TRUE;
509 SCIP_Real* coefs;
510 SCIP_Bool success = FALSE;
511
512 /* if no IPOPT available, don't run test */
513 if( nlpi == NULL )
514 return;
515
516 /* create compicated nlrow */
517 createNlRow3(convex);
518
519 /** point where to compute gradient cut **/
520 SCIP_CALL( SCIPcreateSol(scip, &x0, NULL) );
521 SCIP_CALL( SCIPsetSolVal(scip, x0, x, -0.5) );
522 SCIP_CALL( SCIPsetSolVal(scip, x0, y, 0.5) );
523
524 /* compute gradient cut */
525 SCIP_CALL( SCIPcreateEmptyRowUnspec(scip, &gradcut, "gradcut", -SCIPinfinity(scip), SCIPinfinity(scip),
526 FALSE, FALSE , FALSE) );
527 SCIP_CALL( SCIPexprintCreate(SCIPblkmem(scip), &exprint) );
528 SCIP_CALL( generateCut(scip, x0, exprint, nlrow3, convex, gradcut, &success) );
529 cr_assert(success);
530
531 /* check coefficients */
532 coefs = SCIProwGetVals(gradcut);
533 cr_expect_eq(SCIProwGetNNonz(gradcut), 2);
534 cr_expect_float_eq(-1.29124137035262, coefs[0], EPS, "for x got %f\n", coefs[0]);
535 cr_expect_float_eq(0.307654681677814, coefs[1], EPS, "for y got %f\n", coefs[1]);
536 cr_expect_float_eq(0.936777583155184, SCIProwGetRhs(gradcut), EPS, "wrong rhs\n");
537
538 SCIP_CALL( SCIPexprintFree(&exprint) );
539 SCIP_CALL( SCIPfreeSol(scip, &x0) );
540 SCIP_CALL( SCIPreleaseRow(scip, &gradcut) );
541 }
542
Test(evaluation,gradient_complicated_concave)543 Test(evaluation, gradient_complicated_concave)
544 {
545 SCIP_SOL* x0;
546 SCIP_ROW* gradcut;
547 SCIP_EXPRINT* exprint;
548 SCIP_Bool convex = FALSE;
549 SCIP_Real* coefs;
550 SCIP_Bool success = FALSE;
551
552 /* if no IPOPT available, don't run test */
553 if( nlpi == NULL )
554 return;
555
556 /* create compicated nlrow */
557 createNlRow3(convex);
558
559 /** point where to compute gradient cut **/
560 SCIP_CALL( SCIPcreateSol(scip, &x0, NULL) );
561 SCIP_CALL( SCIPsetSolVal(scip, x0, x, -0.5) );
562 SCIP_CALL( SCIPsetSolVal(scip, x0, y, 0.5) );
563
564 /* compute gradient cut */
565 SCIP_CALL( SCIPcreateEmptyRowUnspec(scip, &gradcut, "gradcut", -SCIPinfinity(scip), SCIPinfinity(scip),
566 FALSE, FALSE , FALSE) );
567 SCIP_CALL( SCIPexprintCreate(SCIPblkmem(scip), &exprint) );
568 SCIP_CALL( generateCut(scip, x0, exprint, nlrow3, convex, gradcut, &success) );
569 cr_assert(success);
570
571 /* check coefficients */
572 coefs = SCIProwGetVals(gradcut);
573 cr_expect_eq(SCIProwGetNNonz(gradcut), 2);
574 cr_expect_float_eq(1.29124137035262, coefs[0], EPS, "for x got %f\n", coefs[0]);
575 cr_expect_float_eq(-0.307654681677814, coefs[1], EPS, "for y got %f\n", coefs[1]);
576 cr_expect_float_eq(-0.936777583155184, SCIProwGetLhs(gradcut), EPS, "wrong rhs\n");
577
578 SCIP_CALL( SCIPexprintFree(&exprint) );
579 SCIP_CALL( SCIPfreeSol(scip, &x0) );
580 SCIP_CALL( SCIPreleaseRow(scip, &gradcut) );
581 }
582
583 /* test that check interior point, notice that it doesn't belong to the test suite evaluation evaluation */
Test(interior_point,compute_interior_point)584 Test(interior_point, compute_interior_point)
585 {
586 SCIP_SEPADATA* sepadata;
587
588 SCIP_CALL( SCIPcreate(&scip) );
589
590 /* include NLPI's */
591 SCIP_CALL( SCIPcreateNlpSolverIpopt(SCIPblkmem(scip), &nlpi) );
592
593 /* if no IPOPT available, don't run test */
594 if( nlpi == NULL )
595 return;
596
597 SCIP_CALL( SCIPincludeNlpi(scip, nlpi) );
598 SCIP_CALL( SCIPincludeExternalCodeInformation(scip, SCIPgetSolverNameIpopt(), SCIPgetSolverDescIpopt()) );
599
600 /* include gauge separator */
601 SCIP_CALL( SCIPincludeSepaGauge(scip) );
602
603 /* create a problem */
604 SCIP_CALL( SCIPcreateProbBasic(scip, "problem") );
605
606 /* change SCIP's stage to be able to create nlrows and rows; ask SCIP to generate an NLP */
607 SCIP_CALL( TESTscipSetStage(scip, SCIP_STAGE_SOLVING, TRUE) );
608
609 /* create and add variables */
610 SCIP_CALL( SCIPcreateVarBasic(scip, &x, "x", -2.5, 3.1, 0.0, SCIP_VARTYPE_CONTINUOUS) );
611 SCIP_CALL( SCIPcreateVarBasic(scip, &y, "y", -2.5, 3.1, 0.0, SCIP_VARTYPE_CONTINUOUS) );
612 SCIP_CALL( SCIPaddVar(scip, x) );
613 SCIP_CALL( SCIPaddVar(scip, y) );
614
615 /* create nlrows, add them to SCIP and release them */
616 createNlRow1(TRUE);
617 createNlRow2(TRUE);
618 createNlRow3(TRUE);
619
620 SCIP_CALL( SCIPaddNlRow(scip, nlrow1) );
621 SCIP_CALL( SCIPaddNlRow(scip, nlrow2) );
622 SCIP_CALL( SCIPaddNlRow(scip, nlrow3) );
623
624 SCIP_CALL( SCIPreleaseNlRow(scip, &nlrow1) );
625 SCIP_CALL( SCIPreleaseNlRow(scip, &nlrow2) );
626 SCIP_CALL( SCIPreleaseNlRow(scip, &nlrow3) );
627
628 /* get sepadata */
629 sepadata = SCIPsepaGetData(SCIPfindSepa(scip, "gauge"));
630
631 /* compute interior point */
632 SCIP_CALL( computeInteriorPoint(scip, sepadata) );
633
634 /* check sepadata stuff changed in call */
635 cr_assert_not(sepadata->skipsepa);
636 cr_assert(sepadata->isintsolavailable);
637 cr_assert_not_null(sepadata->intsol);
638
639 /* check interior solution */
640 cr_expect_float_eq(-0.275971168224138, SCIPgetSolVal(scip, sepadata->intsol, x), 1e-5, "received %.10f instead", SCIPgetSolVal(scip, sepadata->intsol, x));
641 cr_expect_float_eq( 0.318323856389092, SCIPgetSolVal(scip, sepadata->intsol, y), 1e-5, "received %g instead", SCIPgetSolVal(scip, sepadata->intsol, y));
642
643 /* free memory */
644 SCIP_CALL( SCIPreleaseVar(scip, &x) );
645 SCIP_CALL( SCIPreleaseVar(scip, &y) );
646
647 /* free SCIP */
648 SCIP_CALL( SCIPexprintCreate(SCIPblkmem(scip), &sepadata->exprinterpreter) ); /* so that it doesn't complain */
649 SCIP_CALL( SCIPfree(&scip) );
650
651 /* check for memory leaks */
652 cr_assert_eq(BMSgetMemoryUsed(), 0, "There is a memory leak!!");
653 }
654