/home/coin/SVN-release/OS-2.0.0/Couenne/src/branch/infeasibilityVT.cpp

Go to the documentation of this file.
00001 /* $Id: infeasibilityVT.cpp 155 2009-06-16 20:19:39Z pbelotti $
00002  *
00003  * Name:    infeasibilityVT.cpp
00004  * Authors: Pietro Belotti, Carnegie Mellon University
00005  * Purpose: Compute violation transfer of a variable. See Tawarmalani
00006  *          and Sahinidis' book or article on MathProg 2002
00007  *
00008  * (C) Carnegie-Mellon University, 2008.
00009  * This file is licensed under the Common Public License (CPL)
00010  */
00011 
00012 #include "CoinHelperFunctions.hpp"
00013 #include "CouenneProblem.hpp"
00014 #include "CouenneVTObject.hpp"
00015 
00016 
00019 double CouenneVTObject::infeasibility (const OsiBranchingInformation *info, int &way) const {
00020 
00021   int indexVar = reference_ -> Index ();
00022   assert (indexVar >= 0);
00023 
00024   CouNumber tol = CoinMin (COUENNE_EPS, feas_tolerance_);
00025 
00026   // feasible variable
00027   if (info -> upper_ [indexVar] - 
00028       info -> lower_ [indexVar] < tol) {
00029     if (reference_ -> isInteger ()) {
00030       double point = info -> solution_ [reference_ -> Index ()];
00031       if (downEstimate_ <       point  - floor (point)) downEstimate_ =       point  - floor (point);
00032       if (upEstimate_   < ceil (point) -        point)  upEstimate_   = ceil (point) -        point;
00033       return intInfeasibility (point);
00034     } else return (upEstimate_ = downEstimate_ = 0.);
00035   }
00036 
00037   // let Couenne know of current status (variables, bounds)
00038   problem_ -> domain () -> push 
00039     (problem_ -> nVars (),
00040      info -> solution_, 
00041      info -> lower_, 
00042      info -> upper_);
00043 
00044   // debug output ////////////////////////////////////////////////////////////////
00045 
00046   if (jnlst_ -> ProduceOutput (J_DETAILED, J_BRANCHING)) {
00047     printf ("VT infeas on ");
00048     reference_ -> print ();
00049     if (reference_ -> Image ()) { // if no list, print image
00050       printf (" := ");
00051       reference_ -> Image () -> print ();
00052     }
00053     const std::set <int> &dependence = problem_ -> Dependence () [indexVar];
00054     if (dependence.size () > 0) {
00055       printf (" -- ");
00056       for (std::set <int>::const_iterator i = dependence.begin (); 
00057            i != dependence.end (); ++i) {
00058         problem_ -> Var (*i) -> print ();
00059         printf (" ");
00060       }
00061     }
00062     printf ("\n");
00063   }
00064 
00065   // get set of variable indices that depend on reference_
00066   const std::set <int> &dependence = problem_ -> Dependence () [indexVar];
00067 
00068   CouNumber
00069     xcurr    = info -> solution_ [indexVar], // current value of variable
00070     retval   = 0.,                           // will be returned
00071     lFeas    = xcurr,                        // left-most feasible point with all but x_i constant
00072     rFeas    = xcurr,                        // right-most ...
00073     leanLeft = 0.,                           // [0,1],  
00074     fx       = xcurr,                        // value of expression associated with variable (if aux)
00075     maxInf   = 0.;                           // max infeasibility over expressions depending on this
00076 
00077   if (reference_ -> Type () == AUX)
00078     fx = (*(reference_ -> Image ())) ();
00079 
00080   if (dependence.size () == 0) { // this is a top level auxiliary,
00081                                  // nowhere an independent
00082 
00083     // that means, for VT, that letting w vary and keeping everything
00084     // else constant we return the difference between current value
00085     // and value of function at this point
00086 
00087     // these variables may also be disabled in BonCouenneSetup.cpp
00088 
00089     retval = (reference_ -> Type () == AUX) ? // if this is an isolated variable
00090       // check if this w=f(x) is used nowhere and is feasible
00091       (upEstimate_ = downEstimate_ = maxInf = checkInfeasibility (info)) :
00092       0.;
00093 
00094     // check if this w=f(x) is used nowhere and is feasible
00095     retval = upEstimate_ = downEstimate_ = maxInf = checkInfeasibility (info);
00096 
00097   } else {
00098 
00099     // this appears as independent in all auxs of the "dependence" list
00100     // check all non-linear objects containing this variable
00101 
00102     for (std::set <int>::const_iterator i = dependence.begin ();
00103          i != dependence.end (); ++i) {
00104 
00105       const CouenneObject &obj = problem_ -> Objects () [*i];
00106       assert (obj. Reference ());
00107 
00108       CouNumber 
00109         left   = xcurr,
00110         right  = xcurr,
00111         infeas = obj. checkInfeasibility (info);
00112 
00113       if (infeas > maxInf)
00114         maxInf = infeas;
00115 
00116       // check feasibility of nonlinear object
00117       if (infeas > 0.) {
00118 
00119         // measure how far have to go, left and right, to restore
00120         // feasibility (see page 582, Tawarmalani and Sahinidis,
00121         // MathProg A: 99, pp. 563-591)
00122         //if (obj. Reference () -> Image ()) // not needed! obj only has nonlinear objects
00123         obj. Reference () -> Image () -> closestFeasible 
00124           (reference_, obj. Reference (), left, right);
00125 
00126         if (left  < lFeas) lFeas = left;
00127         if (right > rFeas) rFeas = right;
00128       }
00129 
00130       if (jnlst_ -> ProduceOutput (J_MATRIX, J_BRANCHING)) { // debug output
00131         jnlst_ -> Printf (J_MATRIX, J_BRANCHING, "[%g,%g] --> %g - %g = %g (diff = %g - %g = %g): ", 
00132                           left, right, rFeas, lFeas, rFeas - lFeas,
00133                           (obj. Reference () -> Image ()) ? 
00134                           (*(obj. Reference () -> Image ())) () : 0.,  
00135                           (*(obj. Reference ())) (),
00136                           (obj. Reference () -> Image ()) ? 
00137                           (*(obj. Reference () -> Image ())) () - (*(obj. Reference ())) () : 0.);
00138         obj. Reference () -> print (); 
00139         if (obj. Reference () -> Image ()) 
00140           {printf (" := "); obj. Reference() -> Image() -> print();}
00141         printf ("\n");
00142       }
00143     }
00144 
00145     if (lFeas < info -> lower_ [indexVar]) lFeas = info -> lower_ [indexVar];
00146     if (rFeas > info -> upper_ [indexVar]) rFeas = info -> upper_ [indexVar];
00147 
00148     retval = rFeas - lFeas;
00149 
00150     upEstimate_   = rFeas - xcurr;
00151     downEstimate_ = xcurr - lFeas;
00152 
00153     // need a nonzero value for the estimate. Although maxInf is
00154     // related to the expressions that depend on this variable, it is
00155     // still an estimate of how the point will change
00156     if (upEstimate_   <= tol) upEstimate_   = maxInf;
00157     if (downEstimate_ <= tol) downEstimate_ = maxInf;
00158   }
00159 
00161 
00162   // object is infeasible. 
00163   // check how in the middle of the interval we are
00164 
00165   const CouNumber threshold = .5; // can be in [0,0.5]
00166 
00167   if (retval > COUENNE_EPS) {
00168 
00169     leanLeft  = (xcurr - lFeas) / retval;
00170 
00171     if      (leanLeft <     threshold) way = 0;
00172     else if (leanLeft > 1 - threshold) way = 1;
00173     else way = (CoinDrand48 () < 0.5) ? 0 : 1;
00174   }
00175 
00176   // done computing delta. Now transfer violation on LP relaxation ///////////////
00177 
00178   // coefficient in objective is in {0,1} and there is only one
00179   // variable with coefficient 1
00180   CouNumber vt_delta =
00181     (indexVar == problem_ -> Obj (0) -> Body () -> Index ()) ? 1. : 0.;
00182 
00183   for (int i=0, n_el = info -> columnLength_ [indexVar]; i < n_el; i++) {
00184 
00185     int indRow = info -> columnStart_ [indexVar] + i;
00186 
00187     vt_delta += 
00188       fabs (info -> pi_ [info -> row_ [indRow]] * 
00189             info -> elementByColumn_  [indRow]);
00190 
00191     //     vt_delta += 
00192     //       info -> pi_ [info -> row_ [indRow]] * 
00193     //       fabs (info -> elementByColumn_  [indRow]);
00194 
00195     jnlst_ -> Printf (J_MATRIX, J_BRANCHING, "+ (pi[%d]=%g) * (el[%d]=%g) [=%g] --> vtd = %g\n",
00196                       info -> row_ [indRow],
00197                       info -> pi_ [info -> row_ [indRow]], 
00198                       indRow,
00199                       info -> elementByColumn_  [indRow],
00200                       info -> pi_ [info -> row_ [indRow]] * 
00201                       info -> elementByColumn_  [indRow],
00202                       vt_delta);
00203   }
00204 
00205   // weights for VT itself and width of interval
00206   const CouNumber 
00207     alpha = 1.0,
00208     beta  = 0.0;
00209 
00210   jnlst_ -> Printf (J_MATRIX, J_BRANCHING, "return %g * %g + %g * %g + %g * %g --> ", 
00211                     alpha, fabs (retval * vt_delta), beta, retval,
00212                     1-alpha-beta, leanLeft * (1-leanLeft));
00213 
00214   retval = 
00215     alpha          * fabs (retval*vt_delta) +  // violation transfer itself
00216     beta           * retval +                  // width of feasibility interval
00217     (1-alpha-beta) * leanLeft * (1-leanLeft);  // how in the middle of the interval x is
00218 
00219   if (jnlst_ -> ProduceOutput (J_MATRIX, J_BRANCHING)) {
00220     if (retval > CoinMin (COUENNE_EPS, feas_tolerance_)) {
00221       printf ("vt-delta is %-10g [", retval); 
00222       reference_ -> print (); 
00223       if (reference_ -> Image ()) { // if no list, print image
00224         printf (" := ");
00225         reference_ -> Image () -> print ();
00226       }
00227       if (dependence.size () > 0) {
00228         printf (" -- ");
00229         for (std::set <int>::const_iterator i = dependence.begin (); 
00230              i != dependence.end (); ++i) {
00231           problem_ -> Var (*i) -> print ();
00232           printf (" ");
00233         }
00234       }
00235       printf ("]\n");
00236     } else printf ("feasible...\n");
00237   }
00238 
00239   problem_ -> domain () -> pop ();
00240 
00241   if ((retval < tol) &&
00242       (maxInf > tol)) {
00243 
00244     // no improvement with this variable, but need to return nonzero
00245     // if this is an auxiliary whose value is different from relative
00246     // expression
00247     retval = maxInf; 
00248   }
00249 
00250   if (retval < CoinMin (COUENNE_EPS, feas_tolerance_)) 
00251     retval = 0.;
00252 
00253 #define ALMOST_ZERO 1e-8
00254 
00255   assert ((retval < ALMOST_ZERO && upEstimate_ < ALMOST_ZERO && downEstimate_ < ALMOST_ZERO) ||
00256           (retval > ALMOST_ZERO && upEstimate_ > ALMOST_ZERO && downEstimate_ > ALMOST_ZERO));
00257 
00258   return (reference_ -> isInteger ()) ? 
00259     CoinMax (retval, intInfeasibility (info -> solution_ [reference_ -> Index ()])) :
00260     retval;
00261 }

Generated on Mon Aug 3 03:02:19 2009 by  doxygen 1.4.7