00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "CoinHelperFunctions.hpp"
00012
00013 #include "CouenneCutGenerator.hpp"
00014
00015 #include "CouenneExprQuad.hpp"
00016 #include "CouenneProblem.hpp"
00017
00018 using namespace Couenne;
00019
00020
00021
00022 void exprQuad::quadCuts (expression *w, OsiCuts &cs, const CouenneCutGenerator *cg){
00023
00024 #ifdef DEBUG
00025 std::cout<<"Expression has "<< lcoeff_.size () <<" linear terms and "
00026 << nqterms_ <<" quadratic terms." << std::endl;
00027
00028 printf ("Q:");
00029 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00030 int xind = row -> first -> Index ();
00031 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
00032 printf (" <%d,%d,%g>", xind, col -> first -> Index (), col -> second);
00033 }
00034
00035 printf ("\nb:");
00036 for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
00037
00038 printf (" <%d,%g>", el -> first -> Index (), el -> second);
00039
00040 if (c0_)
00041 printf ("; <c0 = %g>", c0_);
00042
00043 printf ("\nBounds: var val lb ub eigval scaled\n");
00044
00045 int index = 0;
00046
00047 for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
00048 i != bounds_.end (); ++i, index++) {
00049
00050 printf ("%3d:\t", index);
00051 i -> first -> print (); printf ("\t");
00052 printf (" %8g [%8g, %8g]",
00053 (*(i -> first)) (), i -> second.first, i -> second.second);
00054
00055 CouNumber
00056 lb = cg -> Problem () -> Lb (i -> first -> Index ()),
00057 ub = cg -> Problem () -> Ub (i -> first -> Index ());
00058
00059 if ((eigen_.size () > 0) &&
00060 (fabs (ub-lb) > COUENNE_EPS))
00061 printf (" --> %8g %8g",
00062 eigen_.begin () -> first,
00063 eigen_.begin () -> first / (ub-lb));
00064
00065 printf ("\n");
00066 }
00067 #endif
00068
00069
00070
00071 CouNumber
00072 varVal = (*w) (),
00073 exprVal = (*this) (),
00074 lambda =
00075 (eigen_.size () == 0) ? 0. :
00076 (varVal < exprVal) ?
00077 CoinMin (0., eigen_.begin () -> first) :
00078 CoinMax (0., eigen_.rbegin () -> first),
00079 convVal = 0.;
00080
00081 enum auxSign sign = cg -> Problem () -> Var (w -> Index ()) -> sign ();
00082
00083
00084 if ((sign == expression::AUX_GEQ && varVal > exprVal) ||
00085 (sign == expression::AUX_LEQ && varVal < exprVal))
00086 return;
00087
00088 const CouenneProblem& problem = *(cg -> Problem ());
00089 const int numcols = problem.nVars ();
00090
00091 const double
00092 *colsol = problem.X (),
00093 *lower = problem.Lb (),
00094 *upper = problem.Ub ();
00095
00096
00097
00098
00099 if (fabs (lambda) > COUENNE_EPS) {
00100
00101 convVal = exprVal;
00102
00103
00104
00105 for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
00106 i != bounds_.end (); ++i) {
00107
00108 int ind = i -> first -> Index ();
00109
00110 CouNumber
00111 xi = colsol [ind],
00112 lb = lower [ind],
00113 ub = upper [ind],
00114 delta = ub-lb;
00115
00116 if (fabs (delta) > COUENNE_EPS)
00117 convVal += lambda * (xi-lb) * (ub-xi) / (delta * delta);
00118 }
00119
00120 if (varVal < exprVal) {if (convVal < varVal) return;}
00121 else {if (convVal > varVal) return;}
00122 }
00123
00124 #ifdef DEBUG
00125 std::cout << "Point to cut: ";
00126 for (int i = 0 ; i < numcols ; i++) std::cout << colsol [i] << ", ";
00127 printf (" (w,f(x),c) = (%g,%g,%g) -- lambda = %g\n", (*w) (), exprVal, convVal, lambda);
00128 #endif
00129
00130
00131 double
00132 *Qxs = new double [numcols],
00133 a0 = -c0_;
00134
00135 CoinFillN (Qxs, numcols, 0.);
00136
00137
00138 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00139
00140 int qi = row -> first -> Index ();
00141
00142 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
00143
00144 int qj = col -> first -> Index ();
00145
00146 CouNumber
00147 qc = col -> second,
00148 xi = colsol [qi],
00149 xj = colsol [qj];
00150
00151 if (qi != qj) {
00152 Qxs [qi] += qc * xj;
00153 Qxs [qj] += qc * xi;
00154 a0 += 2 * qc * xi * xj;
00155 }
00156 else {
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170 a0 += qc * xi * xi;
00171 Qxs [qi] += 2 * qc * xi;
00172 }
00173 }
00174 }
00175
00176 #ifdef DEBUG
00177 printf ("2Qx = ("); for (int i = 0; i < numcols; i++) printf ("%g ", Qxs [i]); printf (")\n");
00178 #endif
00179
00180
00181 for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
00182 Qxs [el -> first -> Index ()] += el -> second;
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 Qxs [w -> Index ()] -= 1;
00196
00197 #ifdef DEBUG
00198 printf ("2Qx = ("); for(int i = 0; i < numcols; i++) printf ("%g ", Qxs [i]); printf (")[%g]\n",a0);
00199 #endif
00200
00201
00202
00203 if (fabs (lambda) > COUENNE_EPS)
00204
00205 for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
00206 i != bounds_.end (); ++i) {
00207
00208 int ind = i -> first -> Index ();
00209
00210 CouNumber
00211 xi = colsol [ind],
00212 lb = lower [ind],
00213 ub = upper [ind],
00214 delta = ub-lb;
00215
00216 if (fabs (delta) > COUENNE_EPS) {
00217
00218 CouNumber normlambda = lambda / (delta*delta),
00219 coeff = normlambda * (lb + ub - 2. * xi);
00220
00221 a0 += normlambda * (lb*ub - xi*xi);
00222
00223
00224
00225 Qxs [ind] += coeff;
00226
00227 }
00228
00229
00230
00231
00232
00233 }
00234
00235
00236 int nnz = 0;
00237 for (int i=0; i < numcols ; i++)
00238 if (fabs (Qxs [i]) > COUENNE_EPS)
00239 nnz++;
00240
00241 #ifdef DEBUG
00242 printf ("2Qx = (");for(int i=0;i<numcols;i++)printf("%g ",Qxs[i]);printf (")[%g], %d nz\n",a0,nnz);
00243 #endif
00244
00245
00246 CoinPackedVector a (false);
00247 a.reserve (nnz);
00248
00249 #ifdef DEBUG
00250 CouNumber lhs = 0, lhsc = 0,
00251 *optimum = cg -> Problem () -> bestSol (),
00252 *current = cg -> Problem () -> X ();
00253 #endif
00254
00255 for (int i=0; i < numcols; i++)
00256
00257 if (fabs (Qxs [i]) > 1.0e-21) {
00258
00259
00260 #ifdef DEBUG
00261 if (optimum) {
00262 printf ("%+g * %g ", Qxs [i], optimum [i]);
00263 lhs += Qxs [i] * optimum [i];
00264 }
00265 lhsc += Qxs [i] * current [i];
00266 #endif
00267 a.insert (i, Qxs [i]);
00268 }
00269
00270 OsiRowCut cut;
00271 cut.setRow (a);
00272
00273 delete [] Qxs;
00274
00275 if (varVal < exprVal) {
00276
00277 cut.setUb (a0);
00278
00279 #ifdef DEBUG
00280 if (optimum && (lhs - a0 > COUENNE_EPS)) {
00281 printf ("cut violates optimal solution: %g > %g\n", lhs, a0);
00282 cut.print ();
00283 }
00284 if (lhsc < a0 + COUENNE_EPS){
00285 printf ("cut (+) is not cutting: ");
00286 cut.print ();
00287 }
00288 #endif
00289
00290 }
00291 else {
00292
00293 cut.setLb (a0);
00294 #ifdef DEBUG
00295 if (optimum && (lhs - a0 < -COUENNE_EPS)) {
00296 printf ("cut violates optimal solution: %g < %g\n", lhs, a0);
00297 cut.print ();
00298 }
00299 if (lhsc > a0 - COUENNE_EPS){
00300 printf ("cut (-) is not cutting: ");
00301 cut.print ();
00302 }
00303 #endif
00304
00305 }
00306
00307 cs.insert (cut);
00308 }