/home/coin/SVN-release/OS-2.4.1/Couenne/src/convex/addEnvelope.cpp

Go to the documentation of this file.
00001 /* $Id: addEnvelope.cpp 560 2011-04-17 10:01:15Z stefan $
00002  *
00003  * File:    addEnvelope.cpp
00004  * Author:  Pietro Belotti, Carnegie Mellon University
00005  * Purpose: Add generic envelope to convex function based on function and its derivative
00006  *
00007  * (C) Carnegie-Mellon University, 2006-10.
00008  * This file is licensed under the Eclipse Public License (EPL)
00009  */
00010 
00011 #include "OsiRowCut.hpp"
00012 #include "CouennePrecisions.hpp"
00013 #include "CouenneTypes.hpp"
00014 #include "CouenneCutGenerator.hpp"
00015 #include "CouenneFunTriplets.hpp"
00016 #include "CouenneProblem.hpp"
00017 #include "CouenneExprVar.hpp"
00018 #include "CoinFinite.hpp"
00019 
00020 using namespace Couenne;
00021 
00022 void CouenneCutGenerator::addEnvelope (OsiCuts &cs, int sign,
00023                                        unary_function f,      // function to be linearized
00024                                        unary_function fprime, // derivative of f
00025                                        int w_ind, int x_ind, 
00026                                        CouNumber x, CouNumber l, CouNumber u,
00027                                        t_chg_bounds *chg,
00028                                        bool is_global) const {
00029 
00030   simpletriplet st (f, fprime, fprime); // don't really care about third parameter...
00031   addEnvelope (cs, sign, &st, w_ind, x_ind, x, l, u, chg, is_global);
00032 }
00033 
00034 
00035 void CouenneCutGenerator::addEnvelope (OsiCuts &cs, int sign,
00036                                        funtriplet *ft,
00037                                        int w_ind, int x_ind, 
00038                                        CouNumber x, CouNumber l, CouNumber u,
00039                                        t_chg_bounds *chg,
00040                                        bool is_global) const {
00041   CouNumber opp_slope = - ft -> Fp (x);
00042 
00043   // TODO: remove check of !firstcall_ if point is available already
00044 
00045   // if bounds are very close, convexify with a single line
00046 
00047   bool cLeft  = !chg || (chg [x_ind].lower() != t_chg_bounds::UNCHANGED) || firstcall_;
00048   bool cRight = !chg || (chg [x_ind].upper() != t_chg_bounds::UNCHANGED) || firstcall_;
00049 
00050   if (fabs (u - l) < COUENNE_EPS) {
00051 
00052     CouNumber x0 = 0.5 * (u+l), fp0 = ft -> Fp (x0);
00053     if (cLeft || cRight) 
00054       createCut (cs, ft -> F(x0) - fp0 * x0, 0, w_ind, 1., x_ind, - fp0);
00055     return;
00056   }
00057 
00058   // Add tangent in any case
00059 
00060   if (((!firstcall_) || ((x >= l) && (x <= u)))
00061       && !CoinIsnan (opp_slope) 
00062       && (fabs (opp_slope) < COUENNE_INFINITY)) {
00063 
00064     if (!(problem_ -> Var (x_ind) -> isInteger ()))
00065 
00066       // normal cut
00067       createCut (cs, ft -> F (x) + opp_slope * x, sign, w_ind, 1., 
00068                  x_ind, opp_slope, -1, 0., is_global);
00069 
00070     else {
00071 
00072       // If the independent variable is integer, the envelope cut can
00073       // be made tighter: it is the line through
00074       //
00075       // (x1,y1) = (floor (x), ft -> F (floor (x))) and
00076       // (x2,y2) = (ceil  (x), ft -> F (ceil  (x)))
00077       //
00078       // therefore the inequality will be
00079       //
00080       // y + ax >=<  b
00081       //
00082       // where the sign is determined by the variable sign, and since
00083       // the line has equation
00084       //
00085       // x  - x1     y  - y1
00086       // -------  =  -------
00087       // x2 - x1     y2 - y1
00088       //
00089       //                  y1 - y2              y1 - y2
00090       // we have that a = ------- and b = y1 + -------  x1
00091       //                  x2 - x1              x2 - x1
00092       //
00093       // Thanks to Sergey for gently encouraging me to do this :-)
00094 
00095       CouNumber x1, x2, y1, y2;
00096 
00097       if ((x1 = floor (x)) < l)
00098         x1 = ceil (l - COUENNE_EPS);
00099 
00100       y1 = ft -> F (x1);
00101 
00102       x2 = ceil (x);
00103 
00104       if (fabs (x2-x1) < COUENNE_EPS)
00105         x2 += 1.;
00106 
00107       y2 = ft -> F (x2);
00108 
00109       if ((x2 > u)   ||
00110           CoinIsnan (y1) || CoinIsnan (y2) || 
00111           !CoinFinite (y1) || !CoinFinite (y2)) // fall back to non-integer cut
00112 
00113         createCut (cs, ft -> F (x) + opp_slope * x, sign, w_ind, 1., 
00114                    x_ind, opp_slope, -1, 0., is_global);
00115 
00116       else {
00117 
00118         CouNumber 
00119           slope = (y1-y2) / (x2-x1),
00120           rhs   = y1 + slope * x1;
00121 
00122         createCut (cs, rhs, sign, w_ind, 1., 
00123                    x_ind, slope, -1, 0., is_global);
00124       }
00125 
00126       // TODO: if the DEPENDENT variable is integer also, in principle
00127       // the cut can be made tighter, but not as simply.
00128 
00129     }
00130   }
00131 
00132   // If this is the first call, add a set of cuts by dividing the
00133   // interval in a grid and add a cut at each point of the grid
00134 
00135   if ((convtype_ == UNIFORM_GRID) || firstcall_) {
00136 
00137     if (cLeft || cRight) {
00138       // now add tangent at each sampling point
00139 
00140       CouNumber 
00141         sample = l, 
00142         step   = (u-l) / nSamples_;
00143 
00144       //    printf ("[%.4f %.4f], step = %.4f, %d samples\n", 
00145       //            l, u, step, nSamples_);
00146 
00147       for (int i = 0; i <= nSamples_; i++) {
00148 
00149         opp_slope = - ft -> Fp (sample);
00150 
00151         if ((fabs (opp_slope) < COUENNE_INFINITY) && 
00152             (fabs (sample-x) > COUENNE_EPS)) { // do not add twice cut at current point
00153 
00154 
00155           // same as above, check integrality of argument
00156 
00157           if (!(problem_ -> Var (x_ind) -> isInteger ()))
00158 
00159             createCut (cs, ft -> F (sample) + opp_slope * sample, sign, 
00160                        w_ind, 1.,
00161                        x_ind, opp_slope, 
00162                        -1, 0., is_global);
00163           else {
00164 
00165             CouNumber x1, x2, y1, y2;
00166 
00167             if ((x1 = floor (sample)) < l)
00168               x1 = ceil (l - COUENNE_EPS);
00169 
00170             y1 = ft -> F (x1);
00171 
00172             x2 = ceil (sample);
00173 
00174             if (fabs (x2-x1) < COUENNE_EPS)
00175               x2 += 1.;
00176 
00177             y2 = ft -> F (x2);
00178 
00179             if ((x2 > u)   ||
00180                 CoinIsnan (y1) || CoinIsnan (y2) || 
00181                 !CoinFinite (y1) || !CoinFinite (y2)) // fall back to non-integer cut
00182 
00183               createCut (cs, ft -> F (sample) + opp_slope * sample, sign, w_ind, 1., 
00184                          x_ind, opp_slope, -1, 0., is_global);
00185           }
00186         }
00187 
00188         //      printf ("  Uniform %d: ", i); cut -> print ();
00189 
00190         sample += step;
00191       }
00192     }
00193   }
00194   else if (convtype_ != CURRENT_ONLY) {
00195 
00196     CouNumber sample = x;
00197 
00198     if (fabs (opp_slope) < COUENNE_INFINITY)
00199       createCut (cs, ft -> F (x) + opp_slope * x, sign, 
00200                  w_ind, 1.,
00201                  x_ind, opp_slope, 
00202                  -1, 0.,
00203                  is_global);
00204       //      printf ("  Current tangent: "); cut -> print ();
00205 
00206     for (int i = 0; i <= nSamples_/2; i++) {
00207 
00208       sample += (x-l) / nSamples_;
00209       opp_slope = - ft -> Fp (sample);
00210 
00211       if (fabs (opp_slope) < COUENNE_INFINITY)
00212         createCut (cs, ft -> F (sample) + opp_slope * sample, sign, 
00213                    w_ind, 1.,
00214                    x_ind, opp_slope, 
00215                    -1, 0.,
00216                    is_global);
00217         //      printf ("  neighbour -%d: ", i); cut -> print ();
00218     }
00219 
00220     sample = x;
00221 
00222     for (int i = 0; i <= nSamples_/2; i++) {
00223 
00224       sample += (u-x) / nSamples_;
00225       opp_slope = - ft -> Fp (sample);
00226 
00227       createCut (cs, ft -> F(sample) + opp_slope * sample, sign, 
00228                  w_ind, 1.,
00229                  x_ind, opp_slope, 
00230                  -1, 0.,
00231                  is_global);
00232         //      printf ("  neighbour  %d: ", i); cut -> print ();
00233     }
00234   }
00235 }

Generated on Thu Nov 10 03:05:43 2011 by  doxygen 1.4.7