/home/coin/SVN-release/OS-2.1.1/Couenne/src/branch/operators/branchExprInv.cpp

Go to the documentation of this file.
00001 /* $Id: branchExprInv.cpp 141 2009-06-03 04:19:19Z pbelotti $ */
00002 /*
00003  * Name:    branchExprInv.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: return branch selection for 1/x
00006  *
00007  * (C) Carnegie-Mellon University, 2006-07.
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include "CoinHelperFunctions.hpp"
00012 
00013 #include "exprInv.hpp"
00014 #include "CouenneObject.hpp"
00015 #include "CouenneBranchingObject.hpp"
00016 #include "projections.hpp"
00017 #include "funtriplets.hpp"
00018 
00019 
00021 CouNumber negPowSelectBranch (const CouenneObject *obj, 
00022                               double * &brpts,  
00023                               double * &brDist, // distance of current LP
00024                                                 // point to new convexifications
00025                               int &way,
00026                               CouNumber k,
00027                               CouNumber x0, CouNumber y0, 
00028                               CouNumber l,  CouNumber u) {
00029 
00030   brDist = (double *) realloc (brDist, 2 * sizeof (double));
00031   brpts  = (double *) realloc (brpts, sizeof (CouNumber));
00032 
00033   // two cases: inside or outside the curves (there are two branches
00034   // of the hyperbola).
00035   //
00036   // Inside: the distance depends on the projection of the current
00037   // point onto the would-be upper envelopes, which forces us to look
00038   // at it numerically. If both bounds are infinite, create a ThreeWay
00039   // branch.
00040   //
00041   // Outside: it suffices to project the current point on the line
00042   // (i.e. to get the closest point on the line) to get the maxi-min
00043   // displacement.
00044   //
00045   // As for all monotonous functions, after choosing *brpts it is
00046   // equivalent to choose w's or x's index as ind, as the implied- and
00047   // propagated bounds will do the rest.
00048 
00049   if ((l < -COUENNE_EPS) && (u > COUENNE_EPS)) { // handle discontinuity
00050 
00051     // no matter if the (negative) exponent is odd or even, we better
00052     // branch on 0 lest we have no good convexification (especially
00053     // with odd exponent)
00054 
00055     *brpts = 0.;
00056     way = TWO_RAND;
00057 
00058     // Closest branch of the hyperbola is on the same side of y+x=0 as
00059     // (x0,y0) => need only one powNewton
00060 
00061     if (fabs (x0) < COUENNE_EPS)
00062       x0 = (x0 <= -0.) ? -COUENNE_EPS : COUENNE_EPS;
00063 
00064     CouNumber xp, xx0 = x0, yy0 = y0, exponent = k;
00065 
00066     // invert dependent and independent if
00067     if ((x0+y0 < 0.) && (x0 > 0.) ||  // in lower half of fourth orthant, or
00068         (x0+y0 > 0.) && (x0 < 0.)) {  // in upper half of second orthant
00069 
00070       exponent = 1. / k;
00071       xx0 = y0;
00072       yy0 = x0;
00073     }
00074 
00075     powertriplet pt (exponent);
00076 
00077     xp = (xx0 >= 0) ? 
00078        powNewton  (xx0,  yy0, &pt) : 
00079       -powNewton (-xx0, -yy0, &pt);
00080 
00081     CouNumber diff = x0 - xp;
00082     y0 -= safe_pow (xp, 1. / k);
00083 
00084     brDist [0] = sqrt (diff*diff + y0*y0); // exact distance
00085     brDist [1] = CoinMax (fabs (x0), 1.);
00086 
00087     if (x0 > 0.) {
00088       double swap = brDist [0];
00089       brDist [0] = brDist [1];
00090       brDist [1] = swap;
00091     }
00092 
00093     return CoinMin (brDist [0], brDist [1]);
00094   }
00095 
00096   int intk = 0;
00097 
00098   bool
00099     isInt    =            fabs (k    - (double) (intk = COUENNE_round (k)))    < COUENNE_EPS,
00100     isInvInt = !isInt && (fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS);
00101 
00102   // case 2: bound interval does not contain zero. Look if inside or
00103   // outside of belly (refer to branchExprExp.cpp)
00104 
00105   if ((x0 >=  0.) &&                          (y0 <  safe_pow  (x0,k))  ||   // x0>0, or
00106       (x0 <= -0.) &&                                                         // x0<0, and
00107       ((isInt &&               !(intk % 2) && (y0 <  safe_pow  (x0,k))) ||     // integer, even
00108        ((isInt || isInvInt) &&  (intk % 2) && (y0 > -safe_pow (-x0,k))))) {    // (inv)integer, odd
00109 
00110     // Inside. Branch on closest point on curve, computed with a
00111     // Newton method
00112 
00113     way = (u < -0.) ? TWO_RIGHT : TWO_LEFT; // explore finite interval first
00114 
00115     powertriplet pt (k);
00116 
00117     *brpts = obj -> midInterval ((x0 >= 0.) ? 
00118                                   powNewton ( x0,  y0, &pt) : 
00119                                  -powNewton (-x0, -y0, &pt), l, u);
00120 
00121     CouNumber dy = y0 - safe_pow (*brpts >= 0 ? *brpts : - *brpts, 1. / k);
00122     x0 -= *brpts;
00123     return (brDist [0] = brDist [1] = sqrt (x0*x0 + dy*dy)); // distance is exact
00124   }
00125 
00126   // Inside, (x0^k) * y0 >= 1. Two cases: /////////////////////////////////////////////////
00127  
00128   // 1) bounds are infinite both horizontally and vertically
00129   // (i.e. [-inf,0] or [0,inf]) --> as for exprExp, pick point on
00130   // diagonal from current to curve, to be sure current will be cut by
00131   // branching rule
00132 
00133   if ((l <   COUENNE_EPS) && (u >   COUENNE_INFINITY) || 
00134       (u > - COUENNE_EPS) && (l < - COUENNE_INFINITY)) {
00135 
00136     /* brpts = (double *) realloc (brpts, 2 * sizeof (double));
00137     way = THREE_CENTER; // focus on central convexification first
00138     brpts [0] = x0;      // draw vertical   from (x0,y0) south (north) to curve y=1/x
00139     brpts [1] = 1. / y0; //      horizontal              west  (east)
00140     CouNumber a = fabs (y0 - 1 / x0), // sides of a triangle with (x0,y0)
00141               b = fabs (x0 - 1 / y0), // as one of the vertices
00142               c = a * cos (atan (a/b)); */
00143 
00144     //brpts = (double *) realloc (brpts, sizeof (double));
00145 
00146     //if (x0 > COUENNE_EPS) 
00147     *brpts = 0.5 * (fabs (x0) + pow (fabs (y0), 1./k));
00148 
00149     if (x0 < 0.) {
00150       *brpts = - *brpts;
00151       brDist [0] = fabs (fabs (y0) - safe_pow (fabs (x0), k));
00152       brDist [1] = *brpts - x0;
00153     } else {
00154       brDist [0] = x0 - *brpts;
00155       brDist [1] = fabs (y0 - safe_pow (x0, k));
00156     }
00157 
00158     //else 
00159     //*brpts = 0.5 * (x0 + pow (y0, 1./k));
00160 
00161     // follow South-East diagonal to find point on curve
00162     // so that current point is surely cut 
00163     //*brpts = 0.5 * (x0 + log (y0)); 
00164     //way = TWO_RAND;
00165     way = (x0 > *brpts) ? TWO_RIGHT : TWO_LEFT;
00166 
00167     return CoinMin (brDist [0], brDist [1]);
00168     //x0 - pow (fabs (y0), 1./k), y0 - pow (x0,k));
00169     //return CoinMin (a, CoinMin (b, c)); // distance is exact
00170   }
00171 
00172   // 2) at least one of them is finite
00173 
00174   if (l < - COUENNE_INFINITY) { // u << -0
00175 
00176     way = TWO_RIGHT;
00177     *brpts = obj -> midInterval (x0, l, u);
00178 
00179     return CoinMin (brDist [0] = y0 - safe_pow (*brpts, 1. / k), 
00180                     brDist [1] = projectSeg (x0, y0, l, safe_pow (l, k), 
00181                                              *brpts, safe_pow (*brpts, k), -1)); // distance is exact
00182   }
00183 
00184   if (u > COUENNE_INFINITY) { // l >> +0
00185 
00186     way = TWO_LEFT;
00187     *brpts = obj -> midInterval (x0, l, u);
00188 
00189     return CoinMin (brDist [1] = y0 - safe_pow (*brpts, 1. / k), 
00190                     brDist [0] = projectSeg (x0, y0, l, safe_pow (l, k), 
00191                                              *brpts, safe_pow (*brpts, k), +1)); // distance is exact
00192   }
00193 
00194   // last case: nice finite interval and limited curve
00195 
00196   powertriplet ft (k);
00197   *brpts = obj -> getBrPoint (&ft, x0, l, u);
00198 
00199   /*  // TODO: check if it works with all exponents
00200   if (u > l + COUENNE_EPS) {
00201 
00202     powertriplet ft (k);
00203     *brpts = maxHeight (&ft, l, u); // min area
00204 
00205     // *brpts = safe_pow ((safe_pow (u,k) - safe_pow (l,k)) / (k * (u-l)), 1/(k-1));
00206     // if (u < 0)
00207     // *brpts = - *brpts;
00208   }
00209   else *brpts = midInterval (x0, l, u);*/
00210 
00211   way = TWO_RAND;
00212 
00213   x0 -= *brpts;
00214   y0 -= safe_pow (*brpts, k);
00215 
00216   brDist [0] = projectSeg (x0,y0, l,      safe_pow (l,      k), *brpts, safe_pow (*brpts, k), 0);
00217   brDist [1] = projectSeg (x0,y0, *brpts, safe_pow (*brpts, k), u,      safe_pow (u,      k), 0);
00218 
00219   return CoinMin (brDist [0], brDist [1]);//sqrt (x0*x0 + y0*y0); // distance is exact
00220 }
00221 
00222 
00223 
00226 
00227 CouNumber exprInv::selectBranch (const CouenneObject *obj, 
00228                                  const OsiBranchingInformation *info,
00229                                  expression *&var,
00230                                  double * &brpts, 
00231                                  double * &brDist, // distance of current LP
00232                                                    // point to new convexifications
00233                                  int &way) {
00234 
00235   var = argument_;
00236 
00237   int
00238     ind = argument_           -> Index (),
00239     wi  = obj -> Reference () -> Index ();
00240 
00241   assert ((ind >= 0) && (wi >= 0));
00242 
00243   CouNumber y0 = info -> solution_ [wi],
00244             x0 = info -> solution_ [ind],
00245             l  = info -> lower_    [ind],
00246             u  = info -> upper_    [ind];
00247 
00248   return negPowSelectBranch (obj, brpts, brDist, way, -1, x0, y0, l,  u);
00249 }

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