/home/coin/SVN-release/OS-2.2.0/Couenne/src/convex/operators/conv-exprSinCos.cpp

Go to the documentation of this file.
00001 /* $Id: conv-exprSinCos.cpp 217 2009-07-08 17:02:07Z pbelotti $
00002  *
00003  * Name:    conv-exprSinCos.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: convexification methods for sines and cosines
00006  *
00007  * (C) Carnegie-Mellon University, 2006-09.
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include <math.h>
00012 #ifndef M_PI
00013 # define M_PI 3.14159265358979323846
00014 #endif
00015 #ifndef M_PI_2
00016 # define M_PI_2 1.57079632679489661923
00017 #endif
00018 
00019 #include "OsiSolverInterface.hpp"
00020 #include "CouenneTypes.hpp"
00021 #include "CouenneCutGenerator.hpp"
00022 #include "exprSin.hpp"
00023 #include "exprCos.hpp"
00024 #include "exprAux.hpp"
00025 
00026 #define NEW_TRIG
00027 
00028 #ifndef NEW_TRIG
00030 int addHexagon (const CouenneCutGenerator *, // cut generator that has called us
00031                 OsiCuts &,                   // cut set to be enriched
00032                 enum cou_trig,               // sine or cosine
00033                 expression *,                // auxiliary variable
00034                 expression *);               // argument of cos/sin (should be a variable)
00035 #endif
00036 
00038 int trigEnvelope (const CouenneCutGenerator *, OsiCuts &,
00039                   expression *, expression *, enum cou_trig);
00040 
00041 
00043 
00044 void exprSin::generateCuts (expression *w, //const OsiSolverInterface &si, 
00045                             OsiCuts &cs, const CouenneCutGenerator *cg,
00046                             t_chg_bounds *chg, int wind, 
00047                             CouNumber lbw, CouNumber ubw) {
00048 
00049   //  int wi = w -> Index ();
00050 
00051   /*if (chg && !(cg -> isFirst ()) && 
00052       (chg [wi].lower() == t_chg_bounds::UNCHANGED) && 
00053       (chg [wi].upper() == t_chg_bounds::UNCHANGED))
00054       return;*/
00055 
00056 #ifdef NEW_TRIG
00057   if (trigEnvelope (cg, cs, w, w -> Image () -> Argument (), COU_SINE) == 0)
00058 #else
00059     if (addHexagon (cg, cs, COU_SINE, w, w -> Image () -> Argument()) == 0)
00060 #endif
00061     {
00062 
00063     }
00064 }
00065 
00066 
00068 
00069 void exprCos::generateCuts (expression *w, //const OsiSolverInterface &si, 
00070                             OsiCuts &cs, const CouenneCutGenerator *cg,
00071                             t_chg_bounds *chg, int wind, 
00072                             CouNumber lbw, CouNumber ubw) {
00073 
00074   //  int wi = w -> Index ();
00075 
00076   /*if (chg && !(cg -> isFirst ()) && 
00077       (chg [wi].lower() == t_chg_bounds::UNCHANGED) && 
00078       (chg [wi].upper() == t_chg_bounds::UNCHANGED))
00079       return;*/
00080 
00081 #ifdef NEW_TRIG
00082   if (trigEnvelope (cg, cs, w, w -> Image () -> Argument (), COU_COSINE) == 0) 
00083 #else
00084     if (addHexagon (cg, cs, COU_COSINE, w, w -> Image () -> Argument()) == 0)
00085 #endif
00086     {
00087 
00088     }
00089 }
00090 
00091 
00093 int bayEnvelope (const CouenneCutGenerator *, OsiCuts &, int, int, 
00094                  CouNumber, CouNumber, CouNumber, bool &, bool &);
00095 
00096 
00098 
00099 int trigEnvelope (const CouenneCutGenerator *cg, // cut generator that has called us
00100                    OsiCuts &cs,                  // cut set to be enriched
00101                    expression *w,
00102                    expression *arg,
00103                    enum cou_trig which_trig) {
00104 
00105   CouNumber lb, ub;
00106   arg -> getBounds (lb, ub);
00107 
00108   // if cosine, scale variables to pretend this is a sine problem
00109   CouNumber displ = (which_trig == COU_COSINE) ? M_PI_2 : 0.;
00110 
00111   int ncuts = 0,
00112     xi = arg -> Index (),
00113     wi = w   -> Index ();
00114 
00115   if (fabs (ub - lb) < COUENNE_EPS) {
00116 
00117     CouNumber x0 = 0.5 * (ub+lb), f, fp;
00118 
00119     if (which_trig == COU_SINE) {f = sin (x0); fp =  cos (x0);}
00120     else                        {f = cos (x0); fp = -sin (x0);}
00121 
00122     return cg -> createCut (cs, f - fp*x0, 0, wi, 1., xi, -fp);
00123   }
00124 
00125   // true if, in the first call (lb), a lower/upper chord was added
00126   // --> no such chord must be generated in the second call (ub)
00127   bool skip_up = false, 
00128        skip_dn = false;
00129 
00130   if (lb > -COUENNE_INFINITY) ncuts += bayEnvelope (cg, cs, wi, xi, lb, ub, displ, skip_up, skip_dn);
00131   if (ub <  COUENNE_INFINITY) ncuts += bayEnvelope (cg, cs, wi, xi, ub, lb, displ, skip_up, skip_dn);
00132 
00133   return ncuts;
00134 }
00135 
00136 
00137 //                             __
00138 // study single bay ( \__/ or /  \ ) of the trigonometric function
00139 //
00140 
00141 int bayEnvelope (const CouenneCutGenerator *cg, // cut generator that has called us
00142                  OsiCuts &cs,                   // cut set to be enriched
00143                  int wi,                        //   dependent variable's index
00144                  int xi,                        // independent
00145                  CouNumber x0,                  // starting point
00146                  CouNumber x1,                  // other bound
00147                  CouNumber displacement,
00148                  bool &skip_up, 
00149                  bool &skip_dn) {
00150 
00151   CouNumber tpt,
00152     rx0  = modulo (x0 + displacement, 2*M_PI),
00153     rx1  = rx0 + x1 - x0,
00154     base = x0 - rx0,
00155     sinrx0 = sin (rx0), zero;
00156 
00157   int ncuts = 0,
00158     up   = (rx0 < M_PI) ? +1 : -1,
00159     left = (x0  < x1)   ? +1 : -1;
00160 
00161   // starting point of the current bay
00162   zero = (up>0) ? 0. : M_PI;
00163 
00164   bool *s0, *s1;
00165 
00166   if (up>0) {s0 = &skip_up; s1 = &skip_dn;}
00167   else      {s0 = &skip_dn; s1 = &skip_up;}
00168 
00169   if (left * (modulo (rx0, M_PI) - M_PI_2) < 0) { 
00170 
00171     // after  flex (i.e., at \_ or /~ ) for left  bound, 
00172     // before flex (i.e., at _/ or ~\ ) for right bound
00173 
00174     // out of the "belly": tangent. If on upper bay consider the lower
00175     // half-plane, and viceversa --> use -up
00176     ncuts += cg -> addTangent (cs, wi, xi, x0, sin (rx0), cos (rx0), -up);
00177 
00178     // leftmost extreme to search for tangent point
00179     CouNumber extr0 = .75 * M_PI * (left+1) - M_PI_2 * up; 
00180 
00181     // in:
00182     if ((left * (rx1 - M_PI * ((left - up) / 2 + 1)) <= 0) ||   // if rx1 in same "belly", or
00183         (left * (rx1 - (tpt = trigNewton
00184                         (rx0, extr0, extr0 + M_PI_2))) <= 0)) { // before closest leaning point 
00185       if (!*s1) // -> chord, if not already added in previous call
00186         *s1 = ((ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1,       sin (rx1), up)) > 0);
00187     } else      ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base+tpt, sin (tpt), up);
00188   }
00189   else {
00190 
00191     // after  stationary point (i.e., _/ or ~\ ) for left  bound, 
00192     // before stationary point (i.e., /~ or \_ ) for right bound
00193   
00194     //    if (left * (rx1 - left * (zero + 5*M_PI_2)) < 0) {
00195     if (left * (rx1 - (4*left - up + 2) * M_PI_2) < 0) {
00196       CouNumber cosrx0 = cos (rx0);
00197       if (up * (sin (rx1) - sinrx0 - cosrx0 * (rx1-rx0)) < 0) 
00198         // (b,sinb) below tangent --> tangent
00199         ncuts += cg -> addTangent (cs, wi, xi, x0, sinrx0, cosrx0, -up);
00200       else {    // up: either chord or leaning plane
00201         CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
00202         tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
00203         if (left * (rx1 - tpt) < 0) {
00204           if (!*s0)
00205             *s0 = ((ncuts += cg->addSegment (cs, wi, xi, x0, sin(rx0), x1,       sin(rx1), -up)) > 0);
00206         } else      ncuts += cg->addSegment (cs, wi, xi, x0, sin(rx0), base+tpt, sin(tpt), -up);
00207       }
00208     } else {
00209       CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
00210       tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
00211       ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), -up);
00212     }
00213 
00214     // down: other chord or leaning plane
00215     if ((left * (rx1 - (zero + M_PI)) < 0) || 
00216         (left * (rx1 - (tpt = trigNewton (rx0, (2 +   left - up) * M_PI_2, 
00217                                                (2 + 2*left - up) * M_PI_2))) < 0)) {
00218       if (!*s1) 
00219         *s1 = ((ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1, sin (rx1), up)) > 0);
00220     } else      ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), up);
00221   }
00222 
00223   return ncuts;
00224 }
00225 
00226 
00227 #ifndef NEW_TRIG
00228 
00229 
00231 
00232 int addHexagon (const CouenneCutGenerator *cg, // cut generator that has called us
00233                 OsiCuts &cs,                   // cut set to be enriched
00234                 enum cou_trig tt,              // sine or cosine
00235                 expression *aux,               // auxiliary variable
00236                 expression *arg) {             // argument of cos/sin (should be a variable)
00237 
00238   // retrieve argument bounds
00239   CouNumber lb, ub;
00240   arg -> getBounds (lb, ub);
00241 
00242   int ncuts = 0,
00243     x_ind = arg -> Index (),
00244     w_ind = aux -> Index ();
00245 
00246   if (fabs (ub - lb) < COUENNE_EPS) {
00247 
00248     CouNumber x0 = 0.5 * (ub+lb), f, fp;
00249 
00250     if (tt == COU_SINE) {f = sin (x0); fp =  cos (x0);}
00251     else                {f = cos (x0); fp = -sin (x0);}
00252 
00253     return cg -> createCut (cs, f - fp*x0, 0, w_ind, 1., x_ind, -fp);
00254   }
00255 
00256   // add  /    \ envelope
00257   //      \    /
00258 
00259   // left
00260   if (lb > -COUENNE_INFINITY) { // if not unbounded
00261     if (tt == COU_SINE) {
00262       ncuts += cg -> createCut (cs, sin (lb) - lb, -1, w_ind, 1., x_ind, -1.); // up: w-x <= f lb - lb
00263       ncuts += cg -> createCut (cs, sin (lb) + lb, +1, w_ind, 1., x_ind,  1.); // dn: w+x >= f lb + lb
00264     }
00265     else {
00266       ncuts += cg -> createCut (cs, cos (lb) - lb, -1, w_ind, 1., x_ind, -1.); // up: w-x <= f lb - lb
00267       ncuts += cg -> createCut (cs, cos (lb) + lb, +1, w_ind, 1., x_ind,  1.); // dn: w+x >= f lb + lb
00268     }
00269   }
00270 
00271   // right
00272   if (ub <  COUENNE_INFINITY) { // if not unbounded
00273     if (tt == COU_SINE) {
00274       ncuts += cg -> createCut (cs, sin (ub) - ub, +1, w_ind, 1., x_ind, -1.); // dn: w-x >= f ub - ub
00275       ncuts += cg -> createCut (cs, sin (ub) + ub, -1, w_ind, 1., x_ind,  1.); // up: w+x <= f ub + ub
00276     }
00277     else {
00278       ncuts += cg -> createCut (cs, cos (ub) - ub, +1, w_ind, 1., x_ind, -1.); // dn: w-x >= f ub - ub
00279       ncuts += cg -> createCut (cs, cos (ub) + ub, -1, w_ind, 1., x_ind,  1.); // up: w+x <= f ub + ub
00280     }
00281   }
00282 
00283   return ncuts;
00284 }
00285 #endif

Generated on Thu Aug 5 03:02:56 2010 by  doxygen 1.4.7