1 /* $Id$
2 *
3 * Name: createCuts.cpp
4 * Author: Pietro Belotti
5 * Purpose: a standard cut creator for use with convexification
6 *
7 * (C) Carnegie-Mellon University, 2006-08.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include "OsiRowCut.hpp"
12
13 #include "CouenneTypes.hpp"
14 #include "CouennePrecisions.hpp"
15 #include "CouenneCutGenerator.hpp"
16 #include "CouenneProblem.hpp"
17
18 using namespace Ipopt;
19 using namespace Couenne;
20
21 /// checks if very large or very small nonzero
badCoeff(CouNumber coe)22 bool badCoeff (CouNumber coe) {
23
24 coe = fabs (coe);
25 return ((coe > COU_MAX_COEFF) || ((coe < COU_MIN_COEFF) && (coe > 0.)));
26 }
27
28
29 /// general procedure for inserting a linear cut with up to three
30 /// variables. Return 1 if cut inserted, 0 if none, <0 if error
31
createCut(OsiCuts & cs,CouNumber lb,CouNumber ub,int i1,CouNumber c1,int i2,CouNumber c2,int i3,CouNumber c3,bool is_global) const32 int CouenneCutGenerator::createCut (OsiCuts &cs,
33 CouNumber lb, CouNumber ub,
34 int i1, CouNumber c1,
35 int i2, CouNumber c2,
36 int i3, CouNumber c3,
37 bool is_global) const {
38 bool numerics = false;
39
40 // a maximum of three terms allowed here. If index -1 (default) the
41 // term is not considered
42
43 int nterms = 0;
44
45 // first check: duplicate index.
46
47 if ((i1 != -1) && (i1 == i3)) {i3 = -1; c1 += c3; c3 = 0;}
48 if ((i2 != -1) && (i2 == i3)) {i3 = -1; c2 += c3; c3 = 0;}
49 if ((i1 != -1) && (i1 == i2)) {i2 = -1; c1 += c2; c2 = 0;}
50
51 // CAUTION: this can make the problem infeasible...
52 if (fabs (c3) <= 1.e-21) { i3 = -1;} // shift coeff/index to
53 if (fabs (c2) <= 1.e-21) { c2 = c3; i2 = i3; i3 = -1;} // keep consistency
54 if (fabs (c1) <= 1.e-21) {c1 = c2; i1 = i2; c2 = c3; i2 = i3; i3 = -1;}
55 // why 1.e-21? Look at CoinPackedMatrix.cpp:2273
56
57 #if 0
58 if (i1 >= 0) {if (fabs (c1) > COU_MAX_COEFF) numerics = true; nterms++;} else c1 = 0;
59 if (i2 >= 0) {if (fabs (c2) > COU_MAX_COEFF) numerics = true; nterms++;} else c2 = 0;
60 if (i3 >= 0) {if (fabs (c3) > COU_MAX_COEFF) numerics = true; nterms++;} else c3 = 0;
61 #else
62 if (i1 >= 0) {if (badCoeff (c1)) numerics = true; nterms++;} else c1 = 0;
63 if (i2 >= 0) {if (badCoeff (c2)) numerics = true; nterms++;} else c2 = 0;
64 if (i3 >= 0) {if (badCoeff (c3)) numerics = true; nterms++;} else c3 = 0;
65 #endif
66
67 if (!nterms) // nonsense cut
68 return 0;
69
70 // cut has large coefficients/rhs, bail out
71 if (numerics
72 //|| ((fabs (lb) < COU_MIN_COEFF) ||
73 //(fabs (ub) < COU_MIN_COEFF))
74 || ((fabs (lb) > COU_MAX_COEFF) &&
75 (fabs (ub) > COU_MAX_COEFF))) {
76
77 jnlst_->Printf(J_STRONGWARNING, J_CONVEXIFYING,
78 "### Discarding cut, large coeff/rhs: %g (%d), %g (%d), %g (%d); [%g,%g]\n",
79 c1, i1, c2, i2, c3, i3, lb, ub);
80 return 0;
81 }
82
83 if (!firstcall_ && addviolated_) { // need to check violation
84
85 const CouNumber *x = problem_ -> X ();
86
87 // compute violation
88 CouNumber violation = 0.;
89
90 if (i1 >= 0) violation += c1 * x [i1];
91 if (i2 >= 0) violation += c2 * x [i2];
92 if (i3 >= 0) violation += c3 * x [i3];
93
94 // quit if not violated
95
96 if ((violation < ub + 0 * COUENNE_EPS) &&
97 (violation > lb - 0 * COUENNE_EPS))
98 return 0;
99 }
100
101 // check if cut violates optimal solution (irrespective of the
102 // branching rules applied, so handle with care)
103
104 CouNumber *best = problem_ -> bestSol ();
105
106 if (best &&
107 ((i1 < 0) || ((best [i1] >= problem_ -> Lb (i1)) && (best [i1] <= problem_ -> Ub (i1)))) &&
108 ((i2 < 0) || ((best [i2] >= problem_ -> Lb (i2)) && (best [i2] <= problem_ -> Ub (i2)))) &&
109 ((i3 < 0) || ((best [i3] >= problem_ -> Lb (i3)) && (best [i3] <= problem_ -> Ub (i3))))) {
110
111 CouNumber lhs = 0.;
112
113 if (i1 >= 0) lhs += c1 * best [i1];
114 if (i2 >= 0) lhs += c2 * best [i2];
115 if (i3 >= 0) lhs += c3 * best [i3];
116
117 if (lhs > ub + COUENNE_EPS)
118 jnlst_->Printf(J_STRONGWARNING, J_CONVEXIFYING,
119 "### cut (%d,%d,%d) (%g,%g,%g) cuts optimum: %g >= %g [%g]\n",
120 i1,i2,i3, c1,c2,c3, lhs, ub, lhs - ub);
121
122 if (lhs < lb - COUENNE_EPS)
123 jnlst_->Printf(J_STRONGWARNING, J_CONVEXIFYING,
124 "### cut (%d,%d,%d) (%g,%g,%g) cuts optimum: %g <= %g [%g]\n",
125 i1,i2,i3, c1,c2,c3, lhs, lb, lb - lhs);
126 }
127
128 // You are here if:
129 //
130 // 1) this is the first call to CouenneCutGenerator::generateCuts(), or
131 // 2) you also want unviolated cuts, or
132 // 3) the cut is violated
133
134 // two cases: cut is of the form w1 [<|>]= alpha, hence a column
135 // cut, or it is of the form (a w1 + b w2 + c w3 [<|>]= alpha), a
136 // row cut
137
138 if ((i2 < 0) && (i3 < 0)) { // column cut /////////////////////////////////////////
139
140 if ( (fabs (c1) < COUENNE_EPS)
141 && (fabs (lb) > COU_MAX_COEFF * COUENNE_EPS)
142 && (fabs (ub) > COU_MAX_COEFF * COUENNE_EPS)) {
143
144 jnlst_->Printf(J_STRONGWARNING, J_CONVEXIFYING,
145 "#### nonsense column cut: %e <= %e w_%d <= %e\n",
146 lb, c1, i1, ub);
147 return 0;
148 }
149
150 OsiColCut *cut = new OsiColCut;
151
152 CouNumber
153 ll = lb / c1,
154 uu = ub / c1;
155
156 if (c1 < 0) {
157 CouNumber tmp = ll;
158 ll = uu;
159 uu = tmp;
160 }
161
162 CouNumber &curL = problem_ -> Lb (i1),
163 &curU = problem_ -> Ub (i1);
164
165 if ((uu < COUENNE_INFINITY) &&
166 (uu < curU - COUENNE_EPS)) {
167
168 cut -> setUbs (1, &i1, &uu);
169 curU = uu; // TODO: chg_bds
170 }
171
172 if ((ll > -COUENNE_INFINITY) &&
173 (ll > curL + COUENNE_EPS)) {
174 cut -> setLbs (1, &i1, &ll);
175 curL = ll; // idem
176 }
177
178 cut -> setGloballyValid (is_global); // global?
179
180 cs.insert (cut);
181 delete cut;
182
183 } else {
184
185 // row cut //////////////////////////////////////////////////////////////////////
186
187 CouNumber *coeff = new CouNumber [nterms];
188 int *index = new int [nterms];
189 OsiRowCut *cut = new OsiRowCut;
190
191 int nt = 0;
192
193 if (i1 >= 0) {coeff [nt] = c1; index [nt++] = i1;}
194 if (i2 >= 0) {coeff [nt] = c2; index [nt++] = i2;}
195 if (i3 >= 0) {coeff [nt] = c3; index [nt++] = i3;}
196
197 if (lb > -COUENNE_INFINITY) cut -> setLb (lb);
198 if (ub < COUENNE_INFINITY) cut -> setUb (ub);
199
200 cut -> setRow (nterms, index, coeff);
201
202 delete [] coeff;
203 delete [] index;
204
205 cut -> setGloballyValid (is_global); // global?
206
207 cs.insert (cut);
208 delete cut;
209 }
210
211 return 1;
212 }
213
214
215 /// general procedure for inserting a linear cut with up to three
216 /// variables. Return 1 if cut inserted, 0 if none, <0 if error
217
createCut(OsiCuts & cs,CouNumber rhs,int sign,int i1,CouNumber c1,int i2,CouNumber c2,int i3,CouNumber c3,bool is_global) const218 int CouenneCutGenerator::createCut (OsiCuts &cs,
219 CouNumber rhs, int sign,
220 int i1, CouNumber c1,
221 int i2, CouNumber c2,
222 int i3, CouNumber c3,
223 bool is_global) const {
224
225 return createCut (cs, (sign >= 0) ? rhs : - COIN_DBL_MAX,
226 (sign <= 0) ? rhs : COIN_DBL_MAX,
227 i1, c1, i2, c2, i3, c3, is_global);
228 }
229