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 convexproj.c
17 * @brief unit test for sepa_convexproj.c methods
18 */
19
20 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
21
22 #include "scip/scip.h"
23 #include "scip/cons_nonlinear.h"
24 #include "nlpi/nlpi_ipopt.h"
25
26 #include "scip/sepa_convexproj.c"
27
28 #include "include/scip_test.h"
29
30 #define EPS 1e-5
31
32 static SCIP* scip = NULL;
33 static SCIP_SEPA* sepa = NULL;
34 static SCIP_NLPI* nlpi = NULL;
35 static SCIP_NLROW* nlrow1 = NULL;
36 static SCIP_NLROW* nlrow2 = NULL;
37 static SCIP_NLROW* nlrow3 = NULL;
38 static SCIP_VAR* x;
39 static SCIP_VAR* y;
40
41 /* This test computes the projection of x = 0.75, y = 0.25 onto
42 * log(exp(x) + exp(y)) <= 1
43 * x^2 <= y
44 * 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
45 * x, y \in [-2.5, 3.1]
46 * and then computes the gradient cut at the projected point.
47 * It does so, by considering each possible combination of representing the constraints,
48 * log(exp(x) + exp(y)) <= 1 or -log(exp(x) + exp(y)) >= -1,
49 * etc.
50 *
51 * In the test we use the following methods:
52 * storeNonlinearConvexNlrows
53 * setQuadraticObj
54 */
55
56 /* create log(exp(x) + exp(y)) <= 1 */
57 static
createNlRow1(SCIP_Bool isrhsconvex)58 void createNlRow1(SCIP_Bool isrhsconvex)
59 {
60 SCIP_EXPRTREE* exprtree;
61 SCIP_EXPRCURV curvature;
62 SCIP_EXPR* xexpr;
63 SCIP_EXPR* yexpr;
64 SCIP_EXPR* exp_x;
65 SCIP_EXPR* exp_y;
66 SCIP_EXPR* sum_exp;
67 SCIP_EXPR* lgexp;
68 SCIP_EXPR* finalexpr;
69 SCIP_Real lhs;
70 SCIP_Real rhs;
71 SCIP_VAR* vars[2];
72
73 /* setup expression for x and y */
74 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &xexpr, SCIP_EXPR_VARIDX, 0) );
75 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &yexpr, SCIP_EXPR_VARIDX, 1) );
76
77 /* setup expression for exp(x) and exp(y) */
78 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exp_x, SCIP_EXPR_EXP, xexpr) );
79 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exp_y, SCIP_EXPR_EXP, yexpr) );
80
81 /* setup expression for exp(x) + exp(y) */
82 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &sum_exp, SCIP_EXPR_PLUS, exp_x, exp_y) );
83
84 /* setup expression for lg(sum) */
85 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &lgexp, SCIP_EXPR_LOG, sum_exp) );
86
87 /* decide curvature */
88 if( isrhsconvex )
89 {
90 curvature = SCIP_EXPRCURV_CONVEX;
91 lhs = -SCIPinfinity(scip);
92 rhs = 1.0;
93 finalexpr = lgexp;
94 }
95 else
96 {
97 SCIP_EXPR* minusone;
98 curvature = SCIP_EXPRCURV_CONCAVE;
99 lhs = -1.0;
100 rhs = SCIPinfinity(scip);
101 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &minusone, SCIP_EXPR_CONST, -1.0) );
102 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &finalexpr, SCIP_EXPR_MUL, minusone, lgexp) );
103 }
104
105 /* setup expression tree with finalexpr as root expression */
106 SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, finalexpr, 2, 0, NULL) );
107
108 /* assign SCIP variables to tree */
109 vars[0] = x;
110 vars[1] = y;
111 SCIP_CALL( SCIPexprtreeSetVars(exprtree, 2, vars) );
112
113 /* create nonlinear row for exprtree <= 1 */
114 SCIP_CALL( SCIPcreateNlRow(scip, &nlrow1, "lg(sum exp)", 0.0, 0, NULL, NULL, 0, NULL, 0, NULL, exprtree, lhs, rhs, curvature) );
115
116 /* release */
117 SCIP_CALL( SCIPexprtreeFree(&exprtree) );
118 }
119
120 /* create x^2 <= y */
121 static
createNlRow2(SCIP_Bool isrhsconvex)122 void createNlRow2(SCIP_Bool isrhsconvex)
123 {
124 SCIP_EXPRTREE* exprtree;
125 SCIP_EXPRCURV curvature;
126 SCIP_EXPR* xexpr;
127 SCIP_EXPR* yexpr;
128 SCIP_EXPR* sqr_x;
129 SCIP_EXPR* diff;
130 SCIP_Real lhs;
131 SCIP_Real rhs;
132 SCIP_VAR* vars[2];
133
134 /* setup expression for x and y */
135 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &xexpr, SCIP_EXPR_VARIDX, 0) );
136 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &yexpr, SCIP_EXPR_VARIDX, 1) );
137
138 /* setup expression for x^2 */
139 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &sqr_x, SCIP_EXPR_SQUARE, xexpr) );
140
141 /* decide curvature */
142 if( isrhsconvex )
143 {
144 curvature = SCIP_EXPRCURV_CONVEX;
145 lhs = -SCIPinfinity(scip);
146 rhs = 0.0;
147 /* setup expression for x^2 - y */
148 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &diff, SCIP_EXPR_MINUS, sqr_x, yexpr) );
149 }
150 else
151 {
152 curvature = SCIP_EXPRCURV_CONCAVE;
153 lhs = 0.0;
154 rhs = SCIPinfinity(scip);
155 /* setup expression for - x^2 + y */
156 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &diff, SCIP_EXPR_MINUS, yexpr, sqr_x) );
157 }
158
159 /* setup expression trees with exprs as root expression */
160 SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, diff, 2, 0, NULL) );
161
162 /* assign SCIP variables to tree */
163 vars[0] = x;
164 vars[1] = y;
165 SCIP_CALL( SCIPexprtreeSetVars(exprtree, 2, vars) );
166
167 /* create nonlinear row for exprtree <= 0 */
168 SCIP_CALL( SCIPcreateNlRow(scip, &nlrow2, "x^2 ? y", 0.0, 0, NULL, NULL, 0, NULL, 0, NULL, exprtree, lhs, rhs, curvature) );
169
170 /* release */
171 SCIP_CALL( SCIPexprtreeFree(&exprtree) );
172 }
173
174 /* 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 */
175 static
createNlRow3(SCIP_Bool isrhsconvex)176 void createNlRow3(SCIP_Bool isrhsconvex)
177 {
178 SCIP_EXPRTREE* exprtree;
179 SCIP_EXPRCURV curvature;
180 SCIP_EXPR* xexpr;
181 SCIP_EXPR* arg_exp_x;
182 SCIP_EXPR* arg_exp_y;
183 SCIP_EXPR* arg_log;
184 SCIP_EXPR* yexpr;
185 SCIP_EXPR* exp_x;
186 SCIP_EXPR* exp_y;
187 SCIP_EXPR* lgexp;
188 SCIP_EXPR* finalexpr;
189 SCIP_Real lhs;
190 SCIP_Real rhs;
191 SCIP_Real side;
192 SCIP_VAR* vars[2];
193
194 /* expression for x and y */
195 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &xexpr, SCIP_EXPR_VARIDX, 0) );
196 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &yexpr, SCIP_EXPR_VARIDX, 1) );
197
198 /* expression for 0.12x+0.1 and 0.1*y */
199 SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &arg_exp_x, 1, &xexpr, (SCIP_Real[]){0.12}, 0.1) );
200 SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &arg_exp_y, 1, &yexpr, (SCIP_Real[]){0.1}, 0.0) );
201
202 /* expression for exp(0.12*x+0.1) and exp(0.1*y) */
203 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exp_x, SCIP_EXPR_EXP, arg_exp_x) );
204 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exp_y, SCIP_EXPR_EXP, arg_exp_y) );
205
206 /* expression for 0.5*exp(0.12*x+0.1) + 2*exp(0.1*y) + 0.7 */
207 SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &arg_log, 2,
208 (SCIP_EXPR*[]){exp_x, exp_y}, (SCIP_Real[]){0.5, 2.0}, 0.7) );
209
210 /* expression for lg(0.5*exp(0.12*x+0.1) + 2*exp(0.1*y) + 0.7) */
211 SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &lgexp, SCIP_EXPR_LOG, arg_log) );
212
213 /* decide curvature */
214 if( isrhsconvex )
215 {
216 curvature = SCIP_EXPRCURV_CONVEX;
217 side = 1.0;
218 lhs = -SCIPinfinity(scip);
219 rhs = 0.5;
220 }
221 else
222 {
223 curvature = SCIP_EXPRCURV_CONCAVE;
224 side = -1.0;
225 lhs = -0.5;
226 rhs = SCIPinfinity(scip);
227 }
228
229 SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &finalexpr, 1,
230 &lgexp, (SCIP_Real[]){side * 0.2}, 0.0) );
231
232 /* build expression tree with log as root expression */
233 SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, finalexpr, 2, 0, NULL) );
234
235 /* assign SCIP variables to tree */
236 vars[0] = x;
237 vars[1] = y;
238 SCIP_CALL( SCIPexprtreeSetVars(exprtree, 2, vars) );
239
240 SCIP_CALL( SCIPcreateNlRow(scip, &nlrow3, "complicated",
241 0.0, 1, &x, (SCIP_Real[]){side*1.1},
242 2, vars,
243 3,
244 (SCIP_QUADELEM[]){
245 {0, 0, side*2.4}, // 2.4 * x^2
246 {0, 1, side*0.01}, // 0.01*x*y
247 {1, 1, side*0.3}
248 },
249 exprtree, lhs, rhs, curvature) );
250
251
252
253 /* release */
254 SCIP_CALL( SCIPexprtreeFree(&exprtree) );
255 }
256
257 /* setup the sepadata */
258 static
setup_sepadata(void)259 void setup_sepadata(void)
260 {
261 SCIP_SEPADATA* sepadata;
262 SCIP_NLROW* nlrows[3] = {nlrow1, nlrow2, nlrow3};
263
264 sepadata = SCIPsepaGetData(sepa);
265 cr_assert(sepadata != NULL);
266
267 SCIP_CALL( storeNonlinearConvexNlrows(scip, sepadata, nlrows, 3) );
268 cr_assert_eq(sepadata->nnlrows, 3, "error: received %d nlrows", sepadata->nnlrows);
269
270 /* initialize some of the sepadata */
271 SCIP_CALL( SCIPexprintCreate(SCIPblkmem(scip), &sepadata->exprinterpreter) );
272 sepadata->nlpinvars = SCIPgetNVars(scip);
273 cr_assert_eq(sepadata->nlpinvars, 2, "error: received %d vars", sepadata->nlpinvars);
274
275 /* create nlpi problem */
276 sepadata->nlpi = SCIPgetNlpis(scip)[0];
277 cr_assert_not_null(sepadata->nlpi);
278
279 SCIP_CALL( SCIPnlpiCreateProblem(sepadata->nlpi, &sepadata->nlpiprob, "convexproj-nlp-unittest") );
280 SCIP_CALL( SCIPhashmapCreate(&sepadata->var2nlpiidx, SCIPblkmem(scip), sepadata->nlpinvars) );
281 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &sepadata->nlpivars, SCIPgetVars(scip), sepadata->nlpinvars) );
282
283 /* I shouldn't care about the cutoff, just assert that the lp solution satisfies the cutoff bound */
284 SCIP_CALL( SCIPcreateNlpiProb(scip, sepadata->nlpi, nlrows, 3,
285 sepadata->nlpiprob, sepadata->var2nlpiidx, NULL, NULL, SCIPgetCutoffbound(scip), FALSE, TRUE) );
286
287 /* set quadratic part of objective function */
288 SCIP_CALL( setQuadraticObj(scip, sepadata) );
289 }
290
291 /* setup of test */
292 static
test_setup(void)293 void test_setup(void)
294 {
295 SCIP_VAR* auxx;
296
297 SCIP_CALL( SCIPcreate(&scip) );
298
299 /* include NLPI's */
300 SCIP_CALL( SCIPcreateNlpSolverIpopt(SCIPblkmem(scip), &nlpi) );
301 /* if no IPOPT available, don't run test */
302 if( nlpi == NULL )
303 return;
304
305 SCIP_CALL( SCIPincludeNlpi(scip, nlpi) );
306 SCIP_CALL( SCIPincludeExternalCodeInformation(scip, SCIPgetSolverNameIpopt(), SCIPgetSolverDescIpopt()) );
307
308 /* include convexproj separator and get it */
309 SCIP_CALL( SCIPincludeSepaConvexproj(scip) );
310 sepa = SCIPfindSepa(scip, "convexproj");
311 cr_assert(sepa != NULL);
312
313 /* create a problem */
314 SCIP_CALL( SCIPcreateProbBasic(scip, "problem") );
315
316 SCIP_CALL( SCIPcreateVarBasic(scip, &auxx, "x", -2.5, 3.1, 1.0, SCIP_VARTYPE_CONTINUOUS) );
317 SCIP_CALL( SCIPaddVar(scip, auxx) );
318 /* change SCIP's stage to be able to create nlrows and rows */
319 SCIP_CALL( TESTscipSetStage(scip, SCIP_STAGE_SOLVING, FALSE) );
320
321 x = SCIPvarGetTransVar(auxx);
322 SCIP_CALL( SCIPreleaseVar(scip, &auxx) );
323 /* create and add variables */
324 SCIP_CALL( SCIPcreateVarBasic(scip, &y, "y", -2.5, 3.1, 0.0, SCIP_VARTYPE_CONTINUOUS) );
325 SCIP_CALL( SCIPaddVar(scip, y) );
326 }
327
328 static
teardown(void)329 void teardown(void)
330 {
331 //SCIP_CALL( SCIPreleaseVar(scip, &x) );
332 SCIP_CALL( SCIPreleaseVar(scip, &y) );
333 if( nlrow1 != NULL )
334 {
335 SCIP_CALL( SCIPreleaseNlRow(scip, &nlrow1) );
336 }
337 if( nlrow2 != NULL )
338 {
339 SCIP_CALL( SCIPreleaseNlRow(scip, &nlrow2) );
340 }
341 if( nlrow3 != NULL )
342 {
343 SCIP_CALL( SCIPreleaseNlRow(scip, &nlrow3) );
344 }
345
346 /* free SCIP */
347 SCIP_CALL( SCIPfree(&scip) );
348
349 /* check for memory leaks */
350 cr_assert_eq(BMSgetMemoryUsed(), 0, "There is a memory leak!!");
351 }
352
353 static
project(SCIP_Bool * isrhsconvex)354 void project(SCIP_Bool* isrhsconvex)
355 {
356 SCIP_SEPADATA* sepadata;
357 SCIP_SOL* toseparate_sol;
358 SCIP_ROW* gradcut;
359 SCIP_Real* coefs;
360 SCIP_Real maxvio;
361 SCIP_RESULT result;
362
363 test_setup();
364 /* if no IPOPT available, don't run test */
365 if( nlpi == NULL )
366 return;
367
368 /* create the nl rows */
369 createNlRow1(isrhsconvex[0]);
370 createNlRow2(isrhsconvex[1]);
371 createNlRow3(isrhsconvex[2]);
372
373 /* setup nlpi and other data */
374 setup_sepadata();
375 sepadata = SCIPsepaGetData(sepa);
376
377 /** point active only on nlrow3 **/
378 /* set point to separate */
379 SCIP_CALL( SCIPcreateSol(scip, &toseparate_sol, NULL) );
380 SCIP_CALL( SCIPsetSolVal(scip, toseparate_sol, x, 0.75) );
381 SCIP_CALL( SCIPsetSolVal(scip, toseparate_sol, y, 0.25) );
382
383 /* compute violations */
384 SCIP_CALL( computeMaxViolation(scip, sepadata, toseparate_sol, &maxvio) );
385 cr_expect_float_eq(0.224076984, sepadata->constraintviolation[0], EPS, "got %f\n", sepadata->constraintviolation[0]);
386 cr_expect_float_eq(0.3125, sepadata->constraintviolation[1], EPS, "got %f\n", sepadata->constraintviolation[1]);
387 cr_expect_float_eq(1.937730557, sepadata->constraintviolation[2], EPS, "got %f\n", sepadata->constraintviolation[2]);
388
389 /* project and compute cuts */
390 SCIP_CALL( separateCuts(scip, sepa, toseparate_sol, &result) );
391
392 /* check cut */
393 cr_assert_eq(result, SCIP_SEPARATED, "result is %d instead of SEPARATED", result);
394 cr_assert_eq(SCIPgetNCuts(scip), 1, "got %d cuts", SCIPgetNCuts(scip));
395 gradcut = SCIPgetCuts(scip)[0];
396 coefs = SCIProwGetVals(gradcut);
397 if( isrhsconvex[2] )
398 {
399 cr_expect_float_eq(1.900159674905818, coefs[0], EPS, "for x got %f\n", coefs[0]);
400 cr_expect_float_eq(0.138455761838885, coefs[1], EPS, "for y got %f\n", coefs[1]);
401 cr_expect_float_eq(0.343035343633182, SCIProwGetRhs(gradcut), EPS, "wrong rhs got %f\n", SCIProwGetRhs(gradcut));
402 }
403 else
404 {
405 cr_expect_float_eq(-1.900159674905818, coefs[0], EPS, "for x got %f\n", coefs[0]);
406 cr_expect_float_eq(-0.138455761838885, coefs[1], EPS, "for y got %f\n", coefs[1]);
407 cr_expect_float_eq(-0.343035343633182, SCIProwGetLhs(gradcut), EPS, "wrong rhs got %f\n", SCIProwGetLhs(gradcut));
408 }
409
410 /* to remove the added cuts */
411 SCIP_CALL( SCIPremoveInefficaciousCuts(scip) );
412 SCIP_CALL( SCIPfreeSol(scip, &toseparate_sol) );
413 teardown();
414 }
415
416 TheoryDataPoints(evaluation, convex_is_minus_concave) =
417 {
418 DataPoints(SCIP_Bool, TRUE, FALSE),
419 DataPoints(SCIP_Bool, TRUE, FALSE),
420 DataPoints(SCIP_Bool, TRUE, FALSE)
421 };
422
423 Theory((SCIP_Bool is1convex, SCIP_Bool is2convex, SCIP_Bool is3convex), evaluation, convex_is_minus_concave)
424 {
425 fprintf(stderr, "calling projection with %d %d %d\n",is1convex, is2convex, is3convex);
426 project((SCIP_Bool[]){is1convex, is2convex, is3convex});
427 }
428