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 cons_quadratic 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/cons_quadratic.c"
27
28 #include "include/scip_test.h"
29
30 #define EPS 1e-6
31
32 static SCIP* scip;
33 static SCIP_NLPI* nlpi;
34
35 /* creates scip, problem, includes nonlinear and quadratic cons handlers, and includes NLP */
36 static
setup(void)37 void setup(void)
38 {
39 SCIP_CONSHDLRDATA* conshdlrdata;
40 SCIP_CONSHDLR* conshdlr;
41 SCIP_CALL( SCIPcreate(&scip) );
42
43 /* include quadratic conshdlr (need to include nonlinear) */
44 SCIP_CALL( SCIPincludeConshdlrNonlinear(scip) );
45 SCIP_CALL( SCIPincludeConshdlrQuadratic(scip) );
46
47 /* activate gauge cuts */
48 conshdlr = SCIPfindConshdlr(scip, CONSHDLR_NAME); /* we are including cons_quadraitc */
49 conshdlrdata = SCIPconshdlrGetData(conshdlr);
50 conshdlrdata->gaugecuts = TRUE;
51
52 /* include NLPI's */
53 SCIP_CALL( SCIPcreateNlpSolverIpopt(SCIPblkmem(scip), &nlpi) );
54 /* if no IPOPT available, don't run test */
55 if( nlpi == NULL )
56 return;
57
58 SCIP_CALL( SCIPincludeNlpi(scip, nlpi) );
59 SCIP_CALL( SCIPincludeExternalCodeInformation(scip, SCIPgetSolverNameIpopt(), SCIPgetSolverDescIpopt()) );
60
61 /* create a problem */
62 SCIP_CALL( SCIPcreateProbBasic(scip, "problem") );
63 }
64
65 Test(separation, gauge, .init = setup,
66 .description = "check computation of gauge function of a convex quadratic set"
67 )
68 {
69 SCIP_VAR* xvar;
70 SCIP_VAR* yvar;
71 SCIP_CONS* cons;
72 SCIP_CONSDATA* consdata;
73 SCIP_CONSHDLR* conshdlr;
74 SCIP_SOL* point;
75 SCIP_Bool success;
76 SCIP_Real val;
77
78 /* if no IPOPT available, don't run test */
79 if( nlpi == NULL )
80 return;
81
82 /* create convex quadratic constraint */
83 {
84 SCIP_VAR* vars[2];
85 SCIP_Real vals[2];
86
87 /* create variables */
88 SCIP_CALL( SCIPcreateVarBasic(scip, &xvar, "x", -SCIPinfinity(scip), SCIPinfinity(scip), 1.0, SCIP_VARTYPE_CONTINUOUS) );
89 SCIP_CALL( SCIPcreateVarBasic(scip, &yvar, "y", -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
90
91 SCIP_CALL( SCIPaddVar(scip, xvar) );
92 SCIP_CALL( SCIPaddVar(scip, yvar) );
93
94 /* create quadratic constraint: x^2 + y^2 <= 1 */
95 vars[0] = xvar;
96 vars[1] = yvar;
97
98 vals[0] = 1.0;
99 vals[1] = 1.0;
100
101 SCIP_CALL( SCIPcreateConsBasicQuadratic(scip, &cons, "test", 0, NULL, NULL, 2, vars, vars, vals, -SCIPinfinity(scip), 1.0) );
102
103 /* tell scip constraint is convex */
104 consdata = SCIPconsGetData(cons);
105 consdata->isconvex = TRUE;
106
107 SCIP_CALL( SCIPaddCons(scip, cons) );
108 }
109
110 /* testing */
111 {
112 SCIP_Real expectedarray[2] = {0.0, 0.0};
113
114 /* compute gauge */
115 conshdlr = SCIPconsGetHdlr(cons);
116 SCIP_CALL( computeGauge(scip, conshdlr, cons) );
117
118 /* this could be a single test probably */
119 cr_assert(consdata->isgaugeavailable, "gauge unavailable, pointless to continue");
120 cr_assert_arr_eq(consdata->interiorpoint, expectedarray, 2, "wrong interior point");
121 cr_assert_float_eq(consdata->interiorpointval, 0.0, EPS, "wrong interior point value");
122 cr_assert_arr_eq(consdata->gaugecoefs, expectedarray, 2, "wrong interior point");
123 cr_assert_float_eq(consdata->gaugeconst, 0.0, EPS, "wrong gauge constant");
124
125 /* create sol where to evaluate the gauge: one could parametrize this
126 * or even create a theory out this, because in our case gauge(X) = X/norm(X)^2 */
127 SCIP_CALL( SCIPcreateSol(scip, &point, NULL) );
128 SCIP_CALL( SCIPsetSolVal(scip, point, xvar, 1.0) );
129 SCIP_CALL( SCIPsetSolVal(scip, point, yvar, 1.0) );
130
131 /* test evaluation */
132 SCIP_CALL( evaluateGauge(scip, conshdlr, cons, point, &val, &success) );
133 cr_assert(success, "unsuccessful evaluation of gauge");
134
135 cr_assert_float_eq(val, 1.41421356237309504880, EPS, "wrong gauge evaluation");
136 }
137
138 SCIP_CALL( SCIPfreeSol(scip, &point) );
139 SCIP_CALL( SCIPreleaseCons(scip, &cons) );
140 SCIP_CALL( SCIPreleaseVar(scip, &xvar) );
141 SCIP_CALL( SCIPreleaseVar(scip, &yvar) );
142
143 /* free SCIP */
144 SCIP_CALL( SCIPfree(&scip) );
145
146 /* check for memory leaks */
147 cr_assert_eq(BMSgetMemoryUsed(), 0, "There is are memory leak!!");
148 }
149