Changeset 548


Ignore:
Timestamp:
May 24, 2007 1:14:05 AM (12 years ago)
Author:
pbelotti
Message:

fixed problem in new convexification cuts for expr{Sin,Cos}. Applying obbt more than once, but obbt disabled due to a bad result in ex2_1_7

Location:
branches/Couenne/Couenne
Files:
11 edited

Legend:

Unmodified
Added
Removed
  • branches/Couenne/Couenne/TODO

    r547 r548  
     1- avoid very large coefficients
    12- Bilinear cuts
    23- trisection
     
    1112- recognition of quadratic forms
    1213- q.f. partitioning into PD + nonconvex
     14- (reverse Polish notation)-based evaluation for efficiency
  • branches/Couenne/Couenne/src/convex/CouenneCutGenerator.cpp

    r543 r548  
    5555
    5656  now = CoinCpuTime ();
    57   //problem_ -> print (std::cout);
    58   //printf ("======================================\n");
     57  //  problem_ -> print (std::cout);
     58  //  printf ("======================================\n");
    5959  problem_ -> standardize ();
    6060
     
    6464  septime_ = now;
    6565
    66   //problem_ -> print (std::cout);
    67   //  printf ("======================================\n");
    68   problem_ -> writeMod ("extended-aw.mod", true);
    69   problem_ -> writeMod ("extended-pb.mod", false);
     66  //  problem_ -> print (std::cout);
     67
     68  //problem_ -> writeMod ("extended-aw.mod", true);
     69  //problem_ -> writeMod ("extended-pb.mod", false);
    7070}
    7171
  • branches/Couenne/Couenne/src/convex/boundTightening.cpp

    r539 r548  
    1212#include <CouenneProblem.h>
    1313#include "BonAuxInfos.hpp"
     14
     15// max # bound tightening iterations
     16#define MAX_BT_ITER 10
    1417
    1518
     
    8083  }
    8184
    82   //////////////////////// Bound tightening ///////////////////////////////////
     85  //////////////////////// Bound propagation and implied bounds ////////////////////
    8386
    8487  int   ntightened = 0,
     
    122125    }
    123126
    124 #define MAX_BT_ITER 10
    125 
    126   } while (ntightened || nbwtightened && (niter++ < MAX_BT_ITER));
    127   // need to check if EITHER procedures gave results, as expression
    128   // structure is not a tree any longer.
     127  } while (((ntightened > 0) || (nbwtightened > 0)) && (niter++ < MAX_BT_ITER));
     128  // continue if EITHER procedures gave (positive) results, as
     129  // expression structure is not a tree.
    129130
    130131  return true;
  • branches/Couenne/Couenne/src/convex/generateCuts.cpp

    r547 r548  
    1717#define LARGE_BOUND 9.999e12
    1818
     19// minimum #bound changed in obbt to generate further cuts
     20#define THRES_NBD_CHANGED 1
     21
     22// depth of the BB tree until which obbt is applied at all nodes
     23#define COU_OBBT_CUTOFF_LEVEL 3
     24
     25// maximum number of obbt iterations
     26#define MAX_OBBT_ITER 4
     27
    1928
    2029// translate changed bound sparse array into a dense one
     
    4049                                        OsiCuts &cs,
    4150                                        const CglTreeInfo info) const {
    42 
    43   //  int nInitCuts = cs.sizeRowCuts ();
    44 
    45   /*printf (":::::::::::::::::::::::: level = %d, pass = %d, intree=%d\n Bounds:\n",
    46     info.level, info.pass, info.inTree);*/
     51  int nInitCuts = cs.sizeRowCuts ();
     52
     53  //  printf (":::::::::::::::::::::::: level = %d, pass = %d, intree=%d\n Bounds:\n",
     54  //  info.level, info.pass, info.inTree);
    4755
    4856  Bonmin::BabInfo * babInfo = dynamic_cast <Bonmin::BabInfo *> (si.getAuxiliaryInfo ());
     
    177185  }
    178186
    179 #define COU_OBBT_CUTOFF_LEVEL 3
    180 
    181187  // change tightened bounds through OsiCuts
    182188  if (nchanged)
     
    185191  //#define USE_OBBT
    186192#ifdef USE_OBBT
    187   if ((info.pass == 0) &&
    188       !firstcall_ &&
     193  if ((!firstcall_ || (info.pass > 0)) &&
    189194      (CoinDrand48 () < (double) COU_OBBT_CUTOFF_LEVEL / (info.level + 1))) {
    190195    // apply OBBT at all levels up to the COU_OBBT_CUTOFF_LEVEL-th,
    191196    // and then with probability inversely proportional to the level
    192197
    193 #define THRES_NBD_CHANGED 1
    194 
    195     int nImprov = obbt (si, cs, chg_bds, babInfo);
     198    int nImprov, nIter = 0;
     199    bool repeat = true;
     200
     201    while (repeat &&
     202           (nIter++ < MAX_OBBT_ITER) &&
     203           ((nImprov = obbt (si, cs, chg_bds, babInfo)) > 0))
     204
     205      if (nImprov >= THRES_NBD_CHANGED) {
     206
     207        /// OBBT has given good results, add convexification with
     208        /// improved bounds
     209
     210        sparse2dense (ncols, chg_bds, changed, nchanged);
     211       
     212        genColCuts (si, cs, nchanged, changed);
     213
     214        int nCurCuts = cs.sizeRowCuts ();
     215        //      genRowCuts (si, cs, nchanged, changed);
     216        repeat = nCurCuts < cs.sizeRowCuts (); // reapply only if new cuts available
     217      }
    196218
    197219    if (nImprov < 0)
    198220      goto end_genCuts;
    199 
    200     if (nImprov >= THRES_NBD_CHANGED) {
    201 
    202       //      if (! (boundTightening (si, cs, chg_bds, babInfo)))
    203       //        goto end_genCuts;
    204 
    205       /// OBBT has given good results, add convexification with
    206       /// improved bounds
    207 
    208       sparse2dense (ncols, chg_bds, changed, nchanged);
    209 
    210       genColCuts (si, cs, nchanged, changed);
    211       genRowCuts (si, cs, nchanged, changed);
    212     }
    213221  }
    214222#endif
     
    228236    ntotalcuts_ = nrootcuts_ = cs.sizeRowCuts ();
    229237  }
    230   else ntotalcuts_ += cs.sizeRowCuts ();
    231                                        
     238  else ntotalcuts_ += (cs.sizeRowCuts () - nInitCuts);
     239
    232240  septime_ += CoinCpuTime () - now;
    233241}
  • branches/Couenne/Couenne/src/convex/obbt.cpp

    r547 r548  
    1212#include <CouenneProblem.h>
    1313
     14#define OBBT_EPS 1e-3
    1415
    1516/// reoptimize and change bound of a variable if needed
    1617bool updateBound (OsiSolverInterface *csi, /// interface to use as a solver
    1718                  int sense,               /// 1: minimize, -1: maximize
    18                   CouNumber &bound) {      /// bound to be updated
     19                  CouNumber &bound,        /// bound to be updated
     20                  bool isint) {            /// is this variable integer
    1921
    20   double opt;
    21 
    22   csi -> setObjSense  (sense); // MINIMIZE
     22  csi -> setObjSense (sense);
    2323  csi -> resolve ();
    2424
    25   if (csi -> isProvenOptimal () &&
    26       ((sense > 0) && ((opt = csi -> getObjValue ()) > bound + COUENNE_EPS)) ||
    27       ((sense < 0) && ((opt = csi -> getObjValue ()) < bound - COUENNE_EPS))) {
     25  if (csi -> isProvenOptimal ()) {
    2826
    29     bound = opt;
     27    double opt = csi -> getObjValue ();
     28
     29    if      ((sense>0) && (opt > bound+OBBT_EPS)) bound = (isint ? ceil (opt) : (opt-OBBT_EPS));
     30    else if ((sense<0) && (opt < bound-OBBT_EPS)) bound = (isint ? floor(opt) : (opt+OBBT_EPS));
     31    else return false;
     32
    3033    return true;
    3134  }
    3235  else return false;
    3336}
     37
    3438
    3539/// Optimality based bound tightening
     
    3842                               char *chg_bds,
    3943                               Bonmin::BabInfo * babInfo) const {
     44  int
     45    nImprov = 0,
     46    ncols   = si.getNumCols (),
     47    objind  = problem_ -> Obj (0) -> Body () -> Index ();
    4048
    41   int nImprov = 0, ncols = si.getNumCols ();
    4249  double *objcoe = (double *) malloc (ncols * sizeof (double));
    4350
     
    5461
    5562  // for all (original+auxiliary) variables x_i,
    56   for (int i=0; i<ncols; i++) {
     63  for (int i=0; i<ncols; i++)
    5764
    58     CouNumber opt;
    59     bool chg = false;
     65    if (i != objind) { // do not improve objective's bounds
    6066
    61     objcoe [i] = 1;
    62     csi -> setObjective (objcoe);
     67      CouNumber opt;
     68      bool chg = false;
     69      int nOrig = problem_ -> nVars ();
     70      bool isInt = (i < nOrig) ?
     71        (problem_ -> Var (i)       -> isInteger ()) :
     72        (problem_ -> Aux (i-nOrig) -> isInteger ());
    6373
    64     // minimize and then maximize x_i on si.
     74      objcoe [i] = 1;
     75      csi -> setObjective (objcoe);
    6576
    66     chg = updateBound (csi,  1, problem_ -> Lb (i));
    67     chg = updateBound (csi, -1, problem_ -> Ub (i)) || chg;
     77      // minimize and then maximize x_i on si.
    6878
    69     if (chg) {
     79      CouNumber l0 = problem_ -> Lb (i);
     80      CouNumber u0 = problem_ -> Ub (i);
    7081
    71       if (!(boundTightening (si, cs, chg_bds, babInfo)))
    72         return -1; // return value to signal infeasibility
     82      if (updateBound (csi,  1, problem_ -> Lb (i), isInt)) {
     83        csi -> setColLower (i, problem_ -> Lb (i));
     84        chg = true;
     85      }
     86     
     87      if (updateBound (csi, -1, problem_ -> Ub (i), isInt)) {
     88        csi -> setColUpper (i, problem_ -> Ub (i));
     89        chg = true;
     90      }
    7391
    74       chg_bds [i] = 1;
    75       nImprov++;
     92      if (chg) {
     93        /*
     94        printf ("%d: [%g,%g] --> [%g,%g]\n", i,
     95                l0, u0,
     96                problem_ -> Lb (i),
     97                problem_ -> Ub (i));
     98        */
     99        if (!(boundTightening (si, cs, chg_bds, babInfo)))
     100          return -1; // tell caller this is infeasible
     101
     102        chg_bds [i] = 1;
     103        nImprov++;
     104      }
     105
     106      // restore null obj fun
     107      objcoe [i] = 0;
    76108    }
    77 
    78     // restore null obj fun
    79     objcoe [i] = 0;
    80   }
    81109
    82110  delete csi;
  • branches/Couenne/Couenne/src/convex/operators/conv-exprPow.cpp

    r523 r548  
    233233        || (u < - COUENNE_EPS)) &&
    234234        (l > - powThres) &&    // and are finite
    235         (u <   powThres))
     235        (u <   powThres) &&
     236        (fabs (l+u) > COUENNE_EPS)) // bounds are not opposite
    236237
    237238      cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l, k), u, safe_pow (u, k), -sign);
  • branches/Couenne/Couenne/src/convex/operators/conv-exprSinCos.cpp

    r534 r548  
    1616#include <exprAux.h>
    1717
    18 //#define NEW_TRIG
     18#define NEW_TRIG
    1919
    2020/// convex cuts for sine or cosine
     
    159159    rx0  = modulo (x0 + displacement, 2*M_PI),
    160160    rx1  = rx0 + x1 - x0,
    161     base = x0 + displacement - rx0,
     161    base = x0 - rx0,
    162162    sinrx0 = sin (rx0), zero;
    163163
    164164  int
    165     up   = (modulo (rx0, 2*M_PI) < M_PI) ? +1 : -1,
    166     left = (x0 < x1)                     ? +1 : -1;
    167 
     165    //    up   = (modulo (rx0, 2*M_PI) < M_PI) ? +1 : -1,
     166    up   = (rx0 < M_PI) ? +1 : -1,
     167    left = (x0  < x1)   ? +1 : -1;
     168
     169  // starting point of the current bay
    168170  zero = (up>0) ? 0. : M_PI;
    169171
  • branches/Couenne/Couenne/src/expression/operators/exprBCos.h

    r543 r548  
    5454inline CouNumber exprLBCos::operator () () {
    5555
    56   register CouNumber l = (*(arglist_ [0])) ();
    57   register CouNumber u = (*(arglist_ [1])) ();
     56  register CouNumber
     57    l = (*(arglist_ [0])) (),
     58    u = (*(arglist_ [1])) ();
    5859
    59   if ((u - l > 2 * M_PI) ||      // 1) interval spans whole cycle
    60       (floor (l/2/M_PI - 0.5) < // 2) there is a 3/2 pi + 2k pi in between
    61        floor (u/2/M_PI - 0.5)))
     60  CouNumber pi2 = 2 * M_PI;
     61 
     62  if ((u - l > pi2) ||       // 1) interval spans whole cycle
     63      (floor (l/pi2 - 0.5) < // 2) there is a pi + 2k pi in between
     64       floor (u/pi2 - 0.5)))
    6265    return -1;
    6366
     
    107110inline CouNumber exprUBCos::operator () () {
    108111
    109   register CouNumber l = (*(arglist_ [0])) ();
    110   register CouNumber u = (*(arglist_ [1])) ();
     112  register CouNumber
     113    l = (*(arglist_ [0])) (),
     114    u = (*(arglist_ [1])) ();
    111115
    112   if ((u - l > 2 * M_PI) || // 1) interval spans whole cycle
    113       (floor (l/2/M_PI) <   // 2) there is a 3/2 pi + 2k pi in between
    114        floor (u/2/M_PI)))
     116  CouNumber pi2 = 2 * M_PI;
     117
     118  if ((u - l > pi2) || // 1) interval spans whole cycle
     119      (floor (l/pi2) < // 2) there is a 3/2 pi + 2k pi in between
     120       floor (u/pi2)))
    115121    return 1;
    116122
  • branches/Couenne/Couenne/src/expression/operators/exprBSin.h

    r543 r548  
    5454inline CouNumber exprLBSin::operator () () {
    5555
    56   register CouNumber l = (*(arglist_ [0])) ();
    57   register CouNumber u = (*(arglist_ [1])) ();
     56  register CouNumber
     57    l = (*(arglist_ [0])) (),
     58    u = (*(arglist_ [1])) ();
    5859
    59   if ((u - l > 2 * M_PI) ||      // 1) interval spans whole cycle
    60       (floor (l/2/M_PI - 0.75) < // 2) there is a 3/2 pi + 2k pi in between
    61        floor (u/2/M_PI - 0.75)))
    62     return -1;
     60  CouNumber pi2 = 2 * M_PI;
     61 
     62  if ((u - l > pi2) ||        // 1) interval spans whole cycle
     63      (floor (l/pi2 - 0.75) < // 2) there is a 3/2 pi + 2k pi in between
     64       floor (u/pi2 - 0.75)))
     65    return -1.;
    6366
    6467  return mymin (sin (l), sin (u));
     
    108111inline CouNumber exprUBSin::operator () () {
    109112
    110   register CouNumber l = (*(arglist_ [0])) ();
    111   register CouNumber u = (*(arglist_ [1])) ();
     113  register CouNumber
     114    l = (*(arglist_ [0])) (),
     115    u = (*(arglist_ [1])) ();
    112116
    113   if ((u - l > 2 * M_PI) ||      // 1) interval spans whole cycle
    114       (floor (l/2/M_PI - 0.25) < // 2) there is a 3/2 pi + 2k pi in between
    115        floor (u/2/M_PI - 0.25)))
    116     return -1;
     117  CouNumber pi2 = 2 * M_PI;
     118
     119  if ((u - l > pi2) ||        // 1) interval spans whole cycle
     120      (floor (l/pi2 - 0.25) < // 2) there is a pi/2 + 2k pi in between
     121       floor (u/pi2 - 0.25)))
     122    return 1.;
    117123
    118124  return mymax (sin (l), sin (u));
  • branches/Couenne/Couenne/src/expression/operators/exprCos.cpp

    r543 r548  
    3434// compute bounds of sin x given bounds of x
    3535void exprCos::getBounds (expression *&lb, expression *&ub) {
    36   lb = new exprConst (-1);
    37   ub = new exprConst (1);
    38   return;
     36  //  lb = new exprConst (-1);
     37  //  ub = new exprConst (1);
     38  //  return;
    3939
    4040  // TODO
  • branches/Couenne/Couenne/src/expression/operators/exprSin.cpp

    r543 r548  
    3838void exprSin::getBounds (expression *&lb, expression *&ub) {
    3939
    40   lb = new exprConst (-1);
    41   ub = new exprConst (1);
    42   return;
     40  //  lb = new exprConst (-1);
     41  //  ub = new exprConst (1);
     42  //  return;
    4343
    4444  // TODO:
     
    4848
    4949  lb = new exprLBSin (xl, xu);
    50   ub = new exprLBSin (new exprClone (xl), new exprClone (xu));
     50  ub = new exprUBSin (new exprClone (xl), new exprClone (xu));
    5151}
Note: See TracChangeset for help on using the changeset viewer.