/home/coin/SVN-release/OS-2.1.1/Couenne/src/bound_tightening/operators/impliedBounds-exprDiv.cpp

Go to the documentation of this file.
00001 /* $Id: impliedBounds-exprDiv.cpp 203 2009-07-07 10:52:11Z pbelotti $ */
00002 /*
00003  * Name:    impliedBounds-exprDiv.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: implied bounds for division operators
00006  *
00007  * (C) Carnegie-Mellon University, 2006. 
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include "exprDiv.hpp"
00012 #include "CouennePrecisions.hpp"
00013 #include "CoinHelperFunctions.hpp"
00014 
00015 
00018 
00019 bool exprDiv::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
00020 
00021   //return false; // !!!
00022 
00023   bool resx, resy = resx = false;
00024 
00025   // y is a constant
00026   if (arglist_ [1] -> Type () == CONST) {
00027 
00028     int ind = arglist_ [0] -> Index ();
00029 
00030     if (ind < 0) {
00031       printf ("exprDiv::impliedBound: Error, w = c/d constants\n");
00032       exit (-1);
00033     }
00034 
00035     CouNumber c = arglist_ [1] -> Value ();
00036 
00037     if (fabs (c) < COUENNE_EPS) {
00038       printf ("exprDiv::impliedBound: Error, division by zero\n");
00039       exit (-1);
00040     }
00041 
00042     // a copy of exprMul::impliedBound for the case where y is a constant
00043 
00044     bool xInt = arglist_ [0] -> isInteger ();
00045 
00046     if (c > COUENNE_EPS) {
00047 
00048       if (updateBound (-1, l+ind, xInt ? ceil (l[wind]*c - COUENNE_EPS) : (l[wind]*c))) {
00049         resx = true; 
00050         chg [ind].setLower (t_chg_bounds::CHANGED);
00051       }
00052 
00053       if (updateBound (+1, u+ind, xInt ? floor (u[wind]*c + COUENNE_EPS) : (u[wind]*c))) {
00054         resx = true; 
00055         chg [ind].setUpper (t_chg_bounds::CHANGED);
00056       }
00057     } 
00058     else if (c < - COUENNE_EPS) {
00059 
00060       if (updateBound (-1, l+ind, xInt ? ceil  (u[wind]*c - COUENNE_EPS) : (u[wind]*c))) {
00061         resx = true; 
00062         chg [ind].setLower (t_chg_bounds::CHANGED);
00063       }
00064 
00065       if (updateBound (+1, u+ind, xInt ? floor (l[wind]*c + COUENNE_EPS) : (l[wind]*c))) {
00066         resx = true; 
00067         chg [ind].setUpper (t_chg_bounds::CHANGED);
00068       }
00069     } 
00070   } else {
00071 
00072     // deal with all other cases
00073 
00074     // Each bound on w is represented on the xy plane with two cones,
00075     // such that joining the extreme rays of both one obtains two
00076     // lines, or in other words, the second cone is obtained through
00077     // the transformation (x',y') = (-x,-y) applied to the first cone.
00078     //
00079     // Bounds can be tightened according to four different cases,
00080     // depending on which corner of the bounding box belongs to which
00081     // of the cones and on the sign of the bounds.
00082     //
00083     // Define wl <= w <= wu,
00084     //        xl <= x <= xu,
00085     //        yl <= y <= yu.
00086     //
00087     // Then the "tightenable" bounds are, depending on the corner:
00088     //
00089     //             _______________________________________________________
00090     //            |           w >= wl         |          w <= wu          |
00091     //            |___________________________|___________________________|
00092     //            |     l<0     |    l>0      |     u<0     |     u>0     |
00093     //            |_____________|_____________|_____________|_____________|
00094     //   Cone     |upper |lower |upper |lower |upper |lower |upper |lower |
00095     //            |      |      |      |      |      |      |      |      |
00096     // 1 xl,yl    | INF  |  -   |  xl  |  yl  |xl,yl?|  yl  |  yl  |  -   |
00097     // 2 xl,yu    |  yu  |  xl  |  -   | INF  |  -   |  yu  |  yu  |yu,xl?|
00098     // 3 xu,yl    |  xu  |  yl  | INF  |  -   |  yl  |  -   |yl,xu?|  yl  |
00099     // 4 xu,yu    |  -   | INF  |  yu  |  xu  |  yu  |xu,yu?|  -   |  yu  |
00100     //            |______|______|______|______|______|______|______|______|
00101     //
00102     // Where "INF" stands for "infeasible subproblem", "-" for
00103     // "nothing to improve", and the rest is improved (those with "?"
00104     // may improve).
00105  
00106     int xi = arglist_ [0] -> Index (),
00107         yi = arglist_ [1] -> Index ();
00108 
00109     CouNumber x0=0;
00110 
00111     CouNumber *xl = l + xi, *yl = l + yi, wl = l [wind],
00112               *xu = u + xi, *yu = u + yi, wu = u [wind];
00113 
00114     /*printf ("from              : w[%d] [%e %e], x%d [%e %e] / y%d [%e %e]",
00115       wind, wl, wu, xi, *xl, *xu, yi, *yl, *yu);*/
00116 
00117     // avoid changing bounds if x is constant
00118     if (xi == -1) 
00119       xl = xu = &x0;
00120 
00122 
00123     // simple case wl = wu = 0
00124 
00125     if ((fabs (wl) < COUENNE_EPS) && 
00126         (fabs (wu) < COUENNE_EPS)) {
00127 
00128       resx = updateBound (-1, xl, 0.) || resx;
00129       resx = updateBound (+1, xl, 0.) || resx;
00130       return resx || resy;
00131     }
00132 
00133     bool resxL, resxU, resyL, 
00134       resyU = resxL = resxU = resyL = false;
00135 
00136     // general case
00137 
00138     if        (wl < - COUENNE_EPS) { // w >= wl, wl negative
00139 
00140       // point C: (xl,yl)
00141 
00142       resyL = ((*yl<0) && (*yl > *xl/wl) && updateBound (-1, yl, 0)) || resyL; //
00143 
00144       if ((*yl>0) && (*yl < *xl/wl)) { // point C violates x/y >= wl, down
00145         resxL = updateBound (-1, xl, *yu*wl) || resxL; //
00146         resyL = updateBound (-1, yl, *xu/wl) || resyL; //
00147       }
00148 
00149       // point B: (xu,yu)
00150 
00151       if ((*yu<0) && (*yu > *xu/wl)) { // point B violates x/y >= wl, down
00152         resxU = updateBound (+1, xu, *yl*wl) || resxU;
00153         resyU = updateBound (+1, yu, *xl/wl) || resyU;
00154       }
00155 
00156       resyU = ((*yu>0) && (*yu < *xu/wl) && updateBound (+1, yu, 0)) || resyU;
00157 
00158     } else if (wl >   COUENNE_EPS) { // w >= wl, wl positive
00159 
00160       //
00161 
00162       resyL = ((*yl<0) && (*yl < *xl/wl) && updateBound (-1, yl, CoinMin (*xl/wl, 0.))) || resyL;
00163       resxL = ((*yl>0) && (*yl > *xl/wl) && updateBound (-1, xl, *yl*wl))               || resxL;
00164 
00165       //
00166 
00167       resxU = ((*yu<0) && (*yu < *xu/wl) && updateBound (+1, xu, *yu*wl))               || resxU;
00168       resyU = ((*yu>0) && (*yu > *xu/wl) && updateBound (+1, yu, CoinMax (*xu/wl, 0.))) || resyU;
00169     }
00170 
00172 
00173     if        (wu >   COUENNE_EPS) { // w <= wu, wu positive
00174 
00175       //
00176 
00177       resyL = ((*yl<0) && (*yl > *xu/wu) && updateBound (-1, yl, 0))      || resyL;
00178 
00179       if ((*yl>0) && (*yl < *xu/wu)) {
00180         resxU = updateBound (+1, xu, *yu*wu) || resxU;
00181         resyL = updateBound (-1, yl, *xl/wu) || resyL;
00182       }
00183 
00184       //
00185 
00186       if ((*yu<0) && (*yu > *xl/wu)) {
00187         resxL = updateBound (-1, xl, *yl*wu) || resxL;
00188         resyU = updateBound (+1, yu, *xu/wu) || resyU;
00189       }
00190 
00191       resyU = ((*yu>0) && (*yu < *xl/wu) && updateBound (+1, yu, 0))      || resyU;
00192 
00193     } else if (wu < - COUENNE_EPS) { // w <= wu, wu negative
00194 
00195       //
00196 
00197       resyL = ((*yl<0) && (*yl < *xu/wu) && updateBound (-1, yl, CoinMin (*xu/wu,0.)))  || resyL;//
00198       resxL = ((*yu<0) && (*yu < *xl/wu) && updateBound (-1, xl, *yu*wu))               || resxL;
00199 
00200       //
00201 
00202       resxU = ((*yl>0) && (*yl > *xu/wu) && updateBound (+1, xu, *yl*wu))               || resxU;
00203       resyU = ((*yu>0) && (*yu > *xl/wu) && updateBound (+1, yu, CoinMax (*xl/wu,0.)))  || resyU;
00204     }
00205 
00206     if (resxL) chg [xi].setLower(t_chg_bounds::CHANGED);
00207     if (resxU) chg [xi].setUpper(t_chg_bounds::CHANGED);
00208     if (resyL) chg [yi].setLower(t_chg_bounds::CHANGED);
00209     if (resyU) chg [yi].setUpper(t_chg_bounds::CHANGED);
00210 
00211     /*if (resx || resy) 
00212       printf ("                 \ntightened division: w[%d] [%e %e], x%d [%e %e] / y%d [%e %e]\n",
00213               wind, wl, wu, xi, *xl, *xu, yi, *yl, *yu);
00214               else printf ("                                                 \r");*/
00215 
00216     resx = resxL || resxU;
00217     resy = resyL || resyU;
00218   }
00219 
00220   bool 
00221     xInt = arglist_ [0] -> isInteger (),
00222     yInt = arglist_ [1] -> isInteger ();
00223 
00224   if (resx && xInt) {
00225     int xi = arglist_ [0] -> Index ();
00226     assert (xi >= 0);
00227     u [xi] = floor (u [xi] + COUENNE_EPS);
00228     l [xi] = ceil  (l [xi] - COUENNE_EPS);
00229   }
00230 
00231   if (resy && yInt) {
00232     int yi = arglist_ [1] -> Index ();
00233     assert (yi >= 0);
00234     u [yi] = floor (u [yi] + COUENNE_EPS);
00235     l [yi] = ceil  (l [yi] - COUENNE_EPS);
00236   }
00237 
00238   return (resx || resy);
00239 }

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