1 /* $Id$
2 *
3 * Name: quadCuts.cpp
4 * Author: Pierre Bonami
5 * Purpose: based on upper and lower convexification, add cuts to convexify
6 *
7 * (C) International Business Machines 2007-10.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include "CoinHelperFunctions.hpp"
12
13 #include "CouenneCutGenerator.hpp"
14
15 #include "CouenneExprQuad.hpp"
16 #include "CouenneProblem.hpp"
17
18 using namespace Couenne;
19
20 //#define DEBUG
21
quadCuts(expression * w,OsiCuts & cs,const CouenneCutGenerator * cg)22 void exprQuad::quadCuts (expression *w, OsiCuts &cs, const CouenneCutGenerator *cg){
23
24 #ifdef DEBUG
25 std::cout<<"Expression has "<< lcoeff_.size () <<" linear terms and "
26 << nqterms_ <<" quadratic terms." << std::endl;
27
28 printf ("Q:");
29 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
30 int xind = row -> first -> Index ();
31 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
32 printf (" <%d,%d,%g>", xind, col -> first -> Index (), col -> second);
33 }
34
35 printf ("\nb:");
36 for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
37 //for (int i=0; i < nlterms_; i++)
38 printf (" <%d,%g>", el -> first -> Index (), el -> second);//index_ [i], coeff_ [i]);
39
40 if (c0_)
41 printf ("; <c0 = %g>", c0_);
42
43 printf ("\nBounds: var val lb ub eigval scaled\n");
44
45 int index = 0;
46
47 for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
48 i != bounds_.end (); ++i, index++) {
49
50 printf ("%3d:\t", index);
51 i -> first -> print (); printf ("\t");
52 printf (" %8g [%8g, %8g]",
53 (*(i -> first)) (), i -> second.first, i -> second.second);
54
55 CouNumber
56 lb = cg -> Problem () -> Lb (i -> first -> Index ()),
57 ub = cg -> Problem () -> Ub (i -> first -> Index ());
58
59 if ((eigen_.size () > 0) &&
60 (fabs (ub-lb) > COUENNE_EPS))
61 printf (" --> %8g %8g",
62 eigen_.begin () -> first,
63 eigen_.begin () -> first / (ub-lb));
64
65 printf ("\n");
66 }
67 #endif
68
69 // Get on which side constraint is violated to get the good lambda
70
71 CouNumber
72 varVal = (*w) (),
73 exprVal = (*this) (),
74 lambda =
75 (eigen_.size () == 0) ? 0. :
76 (varVal < exprVal) ?
77 CoinMin (0., eigen_.begin () -> first) : // Use under-estimator
78 CoinMax (0., eigen_.rbegin () -> first), // Use over-estimator
79 convVal = 0.;
80
81 enum auxSign sign = cg -> Problem () -> Var (w -> Index ()) -> sign ();
82
83 // if this is a "semi"-auxiliary, check if necessary to separate
84 if ((sign == expression::AUX_GEQ && varVal > exprVal) ||
85 (sign == expression::AUX_LEQ && varVal < exprVal))
86 return;
87
88 const CouenneProblem& problem = *(cg -> Problem ());
89 const int numcols = problem.nVars ();
90
91 const double
92 *colsol = problem.X (), // current solution
93 *lower = problem.Lb (), // lower bound
94 *upper = problem.Ub (); // upper
95
96 // compute lower or upper convexification and check if it contains
97 // the current point
98
99 if (fabs (lambda) > COUENNE_EPS) {
100
101 convVal = exprVal;
102
103 // there is a convexification, check if out of current point
104
105 for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
106 i != bounds_.end (); ++i) {
107
108 int ind = i -> first -> Index ();
109
110 CouNumber
111 xi = colsol [ind],
112 lb = lower [ind],
113 ub = upper [ind],
114 delta = ub-lb;
115
116 if (fabs (delta) > COUENNE_EPS)
117 convVal += lambda * (xi-lb) * (ub-xi) / (delta * delta);
118 }
119
120 if (varVal < exprVal) {if (convVal < varVal) return;}
121 else {if (convVal > varVal) return;}
122 }
123
124 #ifdef DEBUG
125 std::cout << "Point to cut: ";
126 for (int i = 0 ; i < numcols ; i++) std::cout << colsol [i] << ", ";
127 printf (" (w,f(x),c) = (%g,%g,%g) -- lambda = %g\n", (*w) (), exprVal, convVal, lambda);
128 #endif
129
130 // Initialize by copying $a$ into a dense vector and computing Q x^*
131 double
132 *Qxs = new double [numcols], // sparse coefficient vector, $Qx^*$
133 a0 = -c0_; // constant term
134
135 CoinFillN (Qxs, numcols, 0.);
136
137 // Compute 2 * Q x^*.
138 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
139
140 int qi = row -> first -> Index ();
141
142 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
143
144 int qj = col -> first -> Index ();
145
146 CouNumber
147 qc = col -> second,
148 xi = colsol [qi],
149 xj = colsol [qj];
150
151 if (qi != qj) {
152 Qxs [qi] += qc * xj; // contribution of element $q_{ij}$ to (Qx)_i
153 Qxs [qj] += qc * xi; // $q_{ij}$ (Qx)_j
154 a0 += 2 * qc * xi * xj;
155 }
156 else {
157 /*
158 if (fabs (lambda) > COUENNE_EPS) {
159
160 CouNumber
161 lb = lower [qi],
162 ub = upper [qi],
163 delta = ub-lb;
164
165 if (fabs (delta) > COUENNE_EPS)
166 qc -= lambda / (delta*delta);
167 }
168 */
169 // elements on the diagonal are not halved upon reading
170 a0 += qc * xi * xi;
171 Qxs [qi] += 2 * qc * xi;
172 }
173 }
174 }
175
176 #ifdef DEBUG
177 printf ("2Qx = ("); for (int i = 0; i < numcols; i++) printf ("%g ", Qxs [i]); printf (")\n");
178 #endif
179
180 // Add a to it.
181 for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
182 Qxs [el -> first -> Index ()] += el -> second; //coeff_ [i];
183
184 // multiply Qx^* by x^*^T again and store the result for the lower
185 // bound into constant term
186
187 /*
188 for (int i=0; i < numcols; i++){
189 a0 -= 0.5 * Qxs [i] * colsol [i];
190 // Qxs [i] *= 2;
191 }
192 */
193
194 // And myself
195 Qxs [w -> Index ()] -= 1;
196
197 #ifdef DEBUG
198 printf ("2Qx = ("); for(int i = 0; i < numcols; i++) printf ("%g ", Qxs [i]); printf (")[%g]\n",a0);
199 #endif
200
201 //a0 -= exprVal;
202
203 if (fabs (lambda) > COUENNE_EPS) // Now the part which depends on lambda, if there is one
204
205 for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
206 i != bounds_.end (); ++i) {
207
208 int ind = i -> first -> Index ();
209
210 CouNumber
211 xi = colsol [ind],
212 lb = lower [ind],
213 ub = upper [ind],
214 delta = ub-lb;
215
216 if (fabs (delta) > COUENNE_EPS) {
217
218 CouNumber normlambda = lambda / (delta*delta),
219 coeff = normlambda * (lb + ub - 2. * xi);
220
221 a0 += normlambda * (lb*ub - xi*xi);
222
223 //a0 += coeff * xi - normlambda * (xi - lb) * (ub - xi);
224 //a0 += normlambda * lb * ub;
225 Qxs [ind] += coeff;
226 //Qxs [ind] += normlambda * (lb + ub);
227 }// else coeff = 0.;
228
229 // a0 += lambda [k] * lower [ind] * upper [ind];
230 // a0 -= lambda [k] * colsol [ind] * colsol [ind];
231
232 //Qxs [ind] -= lambda [k] * (colsol [ind]) * 2;
233 }
234
235 // Count the number of non-zeroes
236 int nnz = 0;
237 for (int i=0; i < numcols ; i++)
238 if (fabs (Qxs [i]) > COUENNE_EPS)
239 nnz++;
240
241 #ifdef DEBUG
242 printf ("2Qx = (");for(int i=0;i<numcols;i++)printf("%g ",Qxs[i]);printf (")[%g], %d nz\n",a0,nnz);
243 #endif
244
245 // Pack the vector into a CoinPackedVector and generate the cut.
246 CoinPackedVector a (false);
247 a.reserve (nnz);
248
249 #ifdef DEBUG
250 CouNumber lhs = 0, lhsc = 0,
251 *optimum = cg -> Problem () -> bestSol (),
252 *current = cg -> Problem () -> X ();
253 #endif
254
255 for (int i=0; i < numcols; i++)
256
257 if (fabs (Qxs [i]) > 1.0e-21) { // why 1.0e-21? Look at CoinPackedMatrix.cpp:2188
258
259 // compute violation
260 #ifdef DEBUG
261 if (optimum) {
262 printf ("%+g * %g ", Qxs [i], optimum [i]);
263 lhs += Qxs [i] * optimum [i];
264 }
265 lhsc += Qxs [i] * current [i];
266 #endif
267 a.insert (i, Qxs [i]);
268 }
269
270 OsiRowCut cut;
271 cut.setRow (a);
272
273 delete [] Qxs;
274
275 if (varVal < exprVal) { //(lambda == dCoeffLo_) {
276
277 cut.setUb (a0);
278
279 #ifdef DEBUG
280 if (optimum && (lhs - a0 > COUENNE_EPS)) {
281 printf ("cut violates optimal solution: %g > %g\n", lhs, a0);
282 cut.print ();
283 }
284 if (lhsc < a0 + COUENNE_EPS){
285 printf ("cut (+) is not cutting: ");
286 cut.print ();
287 }
288 #endif
289 // cut.setLb(-COUENNE_INFINITY);
290 }
291 else {
292
293 cut.setLb (a0);
294 #ifdef DEBUG
295 if (optimum && (lhs - a0 < -COUENNE_EPS)) {
296 printf ("cut violates optimal solution: %g < %g\n", lhs, a0);
297 cut.print ();
298 }
299 if (lhsc > a0 - COUENNE_EPS){
300 printf ("cut (-) is not cutting: ");
301 cut.print ();
302 }
303 #endif
304 // cut.setUb(COUENNE_INFINITY);
305 }
306
307 cs.insert (cut);
308 }
309