/home/coin/SVN-release/OS-2.1.1/Couenne/src/convex/operators/quadCuts.cpp

Go to the documentation of this file.
00001 /* $Id: quadCuts.cpp 217 2009-07-08 17:02:07Z pbelotti $
00002  *
00003  * Name:    quadCuts.cpp
00004  * Author:  Pierre Bonami
00005  * Purpose: based on upper and lower convexification, add cuts to convexify
00006  *
00007  * (C) International Business Machines 2007. 
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include "exprQuad.hpp"
00012 
00013 #include "CouenneProblem.hpp"
00014 #include "CouenneCutGenerator.hpp"
00015 
00016 #include "CoinHelperFunctions.hpp"
00017 
00018 //#define DEBUG
00019 
00020 
00021 void exprQuad::quadCuts (expression *w, OsiCuts &cs, const CouenneCutGenerator *cg){
00022 
00023 #ifdef DEBUG
00024   std::cout<<"Expression has "<< lcoeff_.size () <<" linear terms and "
00025            << nqterms_ <<" quadratic terms." << std::endl;
00026 
00027   printf ("Q:");
00028   for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00029     int xind = row -> first -> Index ();
00030     for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
00031       printf (" <%d,%d,%g>",  xind, col -> first -> Index (), col -> second);
00032   }
00033 
00034   printf ("\nb:");
00035   for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
00036     //for (int i=0; i < nlterms_; i++)
00037     printf (" <%d,%g>",  el -> first -> Index (), el -> second);//index_ [i], coeff_ [i]);
00038 
00039   if (c0_) 
00040     printf ("; <c0 = %g>", c0_);
00041 
00042   printf ("\nBounds: var           val        lb        ub        eigval   scaled\n");
00043 
00044   int index = 0;
00045 
00046   for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
00047        i != bounds_.end (); ++i, index++) {
00048 
00049     printf ("%3d:\t", index);
00050     i -> first -> print (); printf ("\t");
00051     printf (" %8g [%8g, %8g]",
00052             (*(i -> first)) (), i -> second.first, i -> second.second);
00053 
00054     CouNumber 
00055       lb = cg -> Problem () -> Lb (i -> first -> Index ()),
00056       ub = cg -> Problem () -> Ub (i -> first -> Index ());
00057 
00058     if ((eigen_.size () > 0) &&
00059         (fabs (ub-lb) > COUENNE_EPS))
00060       printf (" --> %8g %8g", 
00061               eigen_.begin () -> first,
00062               eigen_.begin () -> first / (ub-lb));
00063 
00064     printf ("\n");
00065   }
00066 #endif
00067 
00068   // Get on which side constraint is violated to get the good lambda
00069 
00070   CouNumber
00071     varVal  = (*w)    (), 
00072     exprVal = (*this) (),
00073     lambda  =
00074     (eigen_.size () == 0) ? 0. :
00075     (varVal < exprVal) ?
00076       CoinMin (0., eigen_.begin  () -> first) : // Use under-estimator
00077       CoinMax (0., eigen_.rbegin () -> first),  // Use  over-estimator
00078     convVal = 0.;
00079 
00080   const CouenneProblem& problem = *(cg -> Problem ());
00081   const int numcols = problem.nVars ();
00082 
00083   const double
00084     *colsol = problem.X  (), // current solution
00085     *lower  = problem.Lb (), //         lower bound
00086     *upper  = problem.Ub (); //         upper
00087 
00088   // compute lower or upper convexification and check if it contains
00089   // the current point
00090 
00091   if (fabs (lambda) > COUENNE_EPS) {
00092 
00093     convVal = exprVal;
00094 
00095     // there is a convexification, check if out of current point
00096 
00097     for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
00098          i != bounds_.end (); ++i) {
00099 
00100       int ind = i -> first -> Index ();
00101 
00102       CouNumber
00103         xi = colsol [ind],
00104         lb = lower  [ind],
00105         ub = upper  [ind],
00106         delta = ub-lb;
00107 
00108       if (fabs (delta) > COUENNE_EPS)
00109         convVal += lambda * (xi-lb) * (ub-xi) / (delta * delta);
00110     }
00111 
00112     if (varVal < exprVal) {if (convVal < varVal) return;}
00113     else                  {if (convVal > varVal) return;}
00114   }
00115 
00116 #ifdef DEBUG
00117   std::cout << "Point to cut: "; 
00118   for (int i = 0 ; i < numcols ; i++) std::cout << colsol [i] << ", ";
00119   printf (" (w,f(x),c) = (%g,%g,%g) -- lambda = %g\n", (*w) (), exprVal, convVal, lambda);
00120 #endif
00121 
00122   // Initialize by copying $a$ into a dense vector and computing Q x^*
00123   double 
00124     *Qxs = new double [numcols], // sparse coefficient vector, $Qx^*$
00125      a0  = -c0_;                 // constant term
00126 
00127   CoinFillN (Qxs, numcols, 0.);
00128 
00129   // Compute 2 * Q x^*.
00130   for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00131 
00132     int qi = row -> first -> Index ();
00133 
00134     for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
00135 
00136       int qj = col -> first -> Index ();
00137 
00138       CouNumber 
00139         qc = col -> second,
00140         xi = colsol [qi],
00141         xj = colsol [qj];
00142 
00143       if (qi != qj) {
00144         Qxs [qi] += qc * xj; // contribution of element $q_{ij}$ to (Qx)_i
00145         Qxs [qj] += qc * xi; //                         $q_{ij}$    (Qx)_j
00146         a0 += 2 * qc * xi * xj;
00147       }
00148       else {
00149         /*
00150         if (fabs (lambda) > COUENNE_EPS) {
00151 
00152           CouNumber
00153             lb = lower  [qi],
00154             ub = upper  [qi],
00155             delta = ub-lb;
00156 
00157           if (fabs (delta) > COUENNE_EPS)
00158             qc -= lambda / (delta*delta);
00159         }
00160         */
00161         // elements on the diagonal are not halved upon reading
00162         a0 += qc * xi * xi;
00163         Qxs [qi] += 2 * qc * xi; 
00164       }
00165     }
00166   }
00167 
00168 #ifdef DEBUG
00169   printf ("2Qx = ("); for (int i = 0; i < numcols; i++) printf ("%g ", Qxs [i]); printf (")\n");
00170 #endif
00171 
00172   // Add a to it.
00173   for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
00174     Qxs [el -> first -> Index ()] += el -> second; //coeff_ [i];
00175 
00176   // multiply Qx^* by x^*^T again and store the result for the lower
00177   // bound into constant term
00178 
00179   /*
00180   for (int i=0; i < numcols; i++){
00181     a0 -= 0.5 * Qxs [i] * colsol [i];
00182     //    Qxs [i] *= 2;
00183   }
00184   */
00185 
00186   // And myself
00187   Qxs [w -> Index ()] -= 1;
00188 
00189 #ifdef DEBUG
00190   printf ("2Qx = ("); for(int i = 0; i < numcols; i++) printf ("%g ", Qxs [i]); printf (")[%g]\n",a0);
00191 #endif
00192 
00193   //a0 -= exprVal;
00194 
00195   if (fabs (lambda) > COUENNE_EPS) // Now the part which depends on lambda, if there is one
00196 
00197     for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
00198          i != bounds_.end (); ++i) {
00199 
00200       int ind = i -> first -> Index ();
00201 
00202       CouNumber 
00203         xi    = colsol [ind],
00204         lb    = lower [ind],
00205         ub    = upper [ind],
00206         delta = ub-lb;
00207 
00208       if (fabs (delta) > COUENNE_EPS) {
00209 
00210         CouNumber normlambda = lambda / (delta*delta),
00211           coeff = normlambda * (lb + ub - 2. * xi);
00212 
00213         a0 += normlambda * (lb*ub - xi*xi);
00214 
00215         //a0 += coeff * xi - normlambda * (xi - lb) * (ub - xi);
00216         //a0 += normlambda * lb * ub;
00217         Qxs [ind] += coeff;
00218         //Qxs [ind] += normlambda * (lb + ub);
00219       }// else coeff = 0.;
00220 
00221       //      a0 += lambda [k] * lower  [ind] * upper  [ind];
00222       //      a0 -= lambda [k] * colsol [ind] * colsol [ind];
00223 
00224       //Qxs [ind] -= lambda [k] * (colsol [ind]) * 2;
00225     }
00226 
00227   // Count the number of non-zeroes
00228   int nnz = 0;
00229   for (int i=0; i < numcols ; i++)
00230     if (fabs (Qxs [i]) > COUENNE_EPS)
00231       nnz++;
00232 
00233 #ifdef DEBUG
00234   printf ("2Qx = (");for(int i=0;i<numcols;i++)printf("%g ",Qxs[i]);printf (")[%g], %d nz\n",a0,nnz);
00235 #endif
00236 
00237   // Pack the vector into a CoinPackedVector and generate the cut.
00238   CoinPackedVector a (false);
00239   a.reserve (nnz);
00240 
00241 #ifdef DEBUG
00242   CouNumber lhs = 0, lhsc = 0,
00243     *optimum = cg -> Problem () -> bestSol (),
00244     *current = cg -> Problem () -> X ();
00245 #endif
00246 
00247   for (int i=0; i < numcols; i++)
00248 
00249     if (fabs (Qxs [i]) > 1.0e-21) { // why 1.0e-21? Look at CoinPackedMatrix.cpp:2188
00250 
00251       // compute violation
00252 #ifdef DEBUG
00253       if (optimum) {
00254         printf ("%+g * %g ", Qxs [i], optimum [i]);
00255         lhs  += Qxs [i] * optimum [i];
00256       }
00257       lhsc += Qxs [i] * current [i];
00258 #endif
00259       a.insert (i, Qxs [i]);
00260     }
00261 
00262   OsiRowCut cut;
00263   cut.setRow (a);
00264 
00265   delete [] Qxs;
00266     
00267   if (varVal < exprVal) { //(lambda == dCoeffLo_) {
00268 
00269      cut.setUb (a0);
00270 
00271 #ifdef DEBUG
00272      if (optimum && (lhs - a0 > COUENNE_EPS)) {
00273        printf ("cut violates optimal solution: %g > %g\n", lhs, a0);
00274        cut.print ();
00275      }
00276      if (lhsc < a0 + COUENNE_EPS){
00277        printf ("cut (+) is not cutting: ");
00278        cut.print ();
00279      }
00280 #endif
00281      //     cut.setLb(-COUENNE_INFINITY);
00282   }
00283   else {
00284 
00285     cut.setLb (a0);
00286 #ifdef DEBUG
00287     if (optimum && (lhs - a0 < -COUENNE_EPS)) {
00288        printf ("cut violates optimal solution: %g < %g\n", lhs, a0);
00289        cut.print ();
00290     }
00291     if (lhsc > a0 - COUENNE_EPS){
00292        printf ("cut (-) is not cutting: ");
00293        cut.print ();
00294      }
00295 #endif
00296     //    cut.setUb(COUENNE_INFINITY);
00297   }
00298 
00299   cs.insert (cut);
00300 }

Generated on Mon May 3 03:05:19 2010 by  doxygen 1.4.7