Changeset 521


Ignore:
Timestamp:
Apr 28, 2007 2:10:25 AM (12 years ago)
Author:
pbelotti
Message:

introducing better bounds for sin/cos. Tool for drawing cuts

Location:
branches/Couenne/Couenne/src
Files:
2 added
10 edited

Legend:

Unmodified
Added
Removed
  • branches/Couenne/Couenne/src/Makefile.am

    r515 r521  
    163163        include/operators/exprPow.h \
    164164        include/operators/exprSin.h \
     165        include/operators/exprBSin.h \
     166        include/operators/exprBCos.h \
    165167        include/operators/exprSub.h \
    166168        include/operators/exprGroup.h \
  • branches/Couenne/Couenne/src/Makefile.in

    r515 r521  
    459459        include/operators/exprPow.h \
    460460        include/operators/exprSin.h \
     461        include/operators/exprBSin.h \
     462        include/operators/exprBCos.h \
    461463        include/operators/exprSub.h \
    462464        include/operators/exprGroup.h \
  • branches/Couenne/Couenne/src/convex/CouenneCutGenerator.cpp

    r497 r521  
    4949  double now = CoinCpuTime ();
    5050
    51   problem_ -> readnl      (asl);
     51  problem_ -> readnl (asl);
    5252
    5353  if ((now = (CoinCpuTime () - now)) > 10.)
  • branches/Couenne/Couenne/src/convex/operators/conv-exprSinCos.cpp

    r497 r521  
    172172  }
    173173  else {
    174     // after stationary point (i.e., _/ or ~\ ) for left bound, before
    175     // for right bound
     174
     175    // after  stationary point (i.e., _/ or ~\ ) for left bound,
     176    // before stationary point (i.e., /~ or \_ ) for right bound
    176177 
    177     if (left * (rx1 - left * (zero + 5*M_PI_2)) < 0) {
     178    //    if (left * (rx1 - left * (zero + 5*M_PI_2)) < 0) {
     179    if (left * (rx1 - (4*left - up + 2) * M_PI_2) < 0) {
    178180      CouNumber cosrx0 = cos (rx0);
    179181      if (up * (sin (rx1) - sinrx0 - cosrx0 * (rx1-rx0)) < 0)
    180182        // (b,sinb) below tangent --> tangent
    181         cg -> addTangent (cs, wi, xi, x0, sinrx0, cosrx0, up);
    182       else     // up: either chord or leaning plane
    183         if (left * (rx1 - (tpt = trigNewton (rx0, left * zero, left * (zero + M_PI_2)))) < 0) {
     183        cg -> addTangent (cs, wi, xi, x0, sinrx0, cosrx0, -up);
     184      else {    // up: either chord or leaning plane
     185        CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
     186        tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
     187        if (left * (rx1 - tpt) < 0) {
    184188          if (!*s0)
    185             *s0 = cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1, sin (rx1), up) > 0;
     189            *s0 = cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1, sin (rx1), -up) > 0;
    186190        }
    187         else cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), up);
    188     } else   cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), up);
     191        else cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), -up);
     192      }
     193    } else {
     194      CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
     195      tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
     196      cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), -up);
     197    }
    189198
    190199    // down: other chord or leaning plane
    191200    if ((left * (rx1 - (zero + M_PI)) < 0) ||
    192         (left * (rx1 - (tpt = trigNewton (rx0, left * zero, left * (zero + M_PI_2)))) < 0)) {
     201        (left * (rx1 - (tpt = trigNewton (rx0,
     202                                          (2 +   left - up) * M_PI_2,
     203                                          (2 + 2*left - up) * M_PI_2))) < 0)) {
    193204      if (!*s1)
    194205        *s1 = (cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1,         sin (rx1), up) > 0);
  • branches/Couenne/Couenne/src/convex/operators/trigNewton.cpp

    r480 r521  
    22 * Name:    trigNewton.cpp
    33 * Author:  Pietro Belotti
    4  * Purpose: numerically find tangents to sines
     4 * Purpose: numerically find tangents to (co)sines
    55 *
    66 * This file is licensed under the Common Public License (CPL)
     
    2525  // F'(x) = - sin x - cos x / (x-a) + (sin x - sin a) / (x - a)^2
    2626
     27  if (l>u) {
     28    register CouNumber swap = l;
     29    l = u;
     30    u = swap;
     31  }
     32
    2733  register CouNumber xk = 0.5 * (u+l);
    2834
     
    4147
    4248    xk += F/Fp;
     49
     50    if      (xk < l) xk = l;
     51    else if (xk > u) xk = u;
    4352
    4453    sinxk = sin (xk);
  • branches/Couenne/Couenne/src/expression/operators/exprCos.cpp

    r331 r521  
    99#include <exprCos.h>
    1010#include <exprSin.h>
     11#include <exprBCos.h>
    1112#include <exprOpp.h>
    1213#include <exprMul.h>
    1314#include <exprClone.h>
     15
    1416#include <math.h>
    1517
    1618
    1719// return an expression -sin (argument), the derivative of cos (argument)
    18 
    1920expression *exprCos::differentiate (int index) {
    2021
     
    2728}
    2829
    29 
    3030// I/O
    31 
    3231void exprCos::print (std::ostream& out) const
    3332  {exprUnary::print (out, "cos", PRE);}
     33
     34// compute bounds of sin x given bounds of x
     35void exprCos::getBounds (expression *&lb, expression *&ub) {
     36  lb = new exprConst (-1);
     37  ub = new exprConst (1);
     38  return;
     39
     40  // TODO
     41  expression *xl, *xu;
     42  argument_ -> getBounds (xl, xu);
     43
     44  lb = new exprLBCos (xl, xu);
     45  ub = new exprUBCos (new exprClone (xl), new exprClone (xu));
     46}
  • branches/Couenne/Couenne/src/expression/operators/exprSin.cpp

    r331 r521  
    11/*
    2  * Name:    exprSin.C
     2 * Name:    exprSin.cpp
    33 * Author:  Pietro Belotti
    44 * Purpose: definition of the sine of a function
     
    1010#include <exprClone.h>
    1111#include <exprCos.h>
     12#include <exprBSin.h>
    1213#include <exprMul.h>
    1314
     
    3132  exprUnary::print (out, "sin", PRE);
    3233}
     34
     35
     36// compute bounds of sin x given bounds of x
     37
     38void exprSin::getBounds (expression *&lb, expression *&ub) {
     39
     40  lb = new exprConst (-1);
     41  ub = new exprConst (1);
     42  return;
     43
     44  // TODO:
     45  expression *xl, *xu;
     46
     47  argument_ -> getBounds (xl, xu);
     48
     49  lb = new exprLBSin (xl, xu);
     50  ub = new exprLBSin (new exprClone (xl), new exprClone (xu));
     51}
  • branches/Couenne/Couenne/src/include/operators/exprCos.h

    r480 r521  
    4141
    4242  /// Get lower and upper bound of an expression (if any)
    43   virtual inline void getBounds (expression *&lb, expression *&ub)
    44     {lb = new exprConst (-1); ub = new exprConst (1);}
     43  void getBounds (expression *&, expression *&);
    4544
    4645  /// generate equality between *this and *w
  • branches/Couenne/Couenne/src/include/operators/exprSin.h

    r477 r521  
    4444
    4545  /// Get lower and upper bound of an expression (if any)
    46   /// TODO: improve with domain
    47   virtual inline void getBounds (expression *&lb, expression *&ub)
    48     {lb = new exprConst (-1); ub = new exprConst (1);}
     46  void getBounds (expression *&, expression *&);
    4947
    5048  /// generate equality between *this and *w
  • branches/Couenne/Couenne/src/util/drawCuts.cpp

    r515 r521  
    2222                   minY =  COUENNE_INFINITY;
    2323
    24   if (1 || (img -> code () == COU_EXPRSIN) ||
     24  if ((img -> code () == COU_EXPRSIN) ||
    2525      (img -> code () == COU_EXPRPOW) ||
    2626      (img -> code () == COU_EXPREXP) ||
     
    2828      (img -> code () == COU_EXPRCOS)) {
    2929
    30     printf (" ==> "); w -> print (std::cout); printf ("\n");
     30    //    fprintf (stderr, " ==> "); w -> print (std::cerr); fprintf (stderr, "\n");
    3131
    3232    expression *lbe, *ube;
     
    3838
    3939    int xi = indep -> Index ();
    40     printf ("looking into w_%d = f (x_%d)\n", w -> Index (), xi);
     40
     41    //    fprintf (stderr, "looking into w_%d = f (x_%d)\n",
     42    //       w -> Index (), xi);
    4143
    4244    indep -> getBounds (lbe, ube);
     
    7173          if (y < minY) minY = y;
    7274
    73           printf ("#=# %.12e %.12e\n", x, y);
     75          fprintf (stderr, "%.12e %.12e\n", x, y);
    7476        }
     77
     78        maxY += (maxY-minY) / 20;
     79        minY -= (maxY-minY) / 20;
    7580      }
    7681       
    77       //lb -= 1;
    78       //ub += 1;
     82      lb -= (ub-lb) / 20;
     83      ub += (ub-lb) / 20;
    7984
    8085      // plot lines defining constraint (only for cuts involving at
     
    9398        }
    9499
    95         printf ("#=# #m=2,S=%d\n", (cs.rowCutPtr (jj) -> sense () == 'L') ? 10:11);
     100        fprintf (stderr, "#m=2,S=%d\n", (cs.rowCutPtr (jj) -> sense () == 'L') ? 10:11);
    96101
    97         printf ("#=# %.12e %.12e\n", lb0, (rhs - el [1] * lb0) / el [0]);
    98         printf ("#=# %.12e %.12e\n", ub0, (rhs - el [1] * ub0) / el [0]);
     102        fprintf (stderr, "%.12e %.12e\n", lb0, (rhs - el [1] * lb0) / el [0]);
     103        fprintf (stderr, "%.12e %.12e\n", ub0, (rhs - el [1] * ub0) / el [0]);
    99104      }
    100105
    101106      cg -> Problem () -> X () [xi] = curx;
     107      exit(0);
    102108    }
    103109  }
Note: See TracChangeset for help on using the changeset viewer.