Changeset 593


Ignore:
Timestamp:
Jun 12, 2007 12:22:32 PM (12 years ago)
Author:
pbelotti
Message:

added warmstart and other checks to OBBT upon suggestion by J. Forrest. Fixed a segfault due to exprMul

Location:
branches/Couenne/Couenne/src/convex
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/Couenne/Couenne/src/convex/obbt.cpp

    r586 r593  
    2424
    2525  csi -> setDblParam (OsiDualObjectiveLimit, COIN_DBL_MAX);
    26   //  csi -> setDblParam (OsiPrimalObjectiveLimit, bound);
    2726  csi -> setDblParam (OsiPrimalObjectiveLimit, (sense==1) ? bound : -bound);
    28 
    29   /*char fname [20];
    30   static int numsave = 0;
    31   sprintf (fname, "obbt_%s_%d", (sense>0) ? "max" : "min", numsave++);
    32   csi -> writeMps (fname);
    33   printf ("\r%s", fname); fflush (stdout);*/
    34 
    3527  csi -> setObjSense (1);
    36   //csi -> initialSolve ();
     28
    3729  csi -> resolve ();
    3830
     
    4032
    4133    double opt = csi -> getObjValue ();
    42     //printf ("found optimum: %g\n", opt);
     34    //    printf ("found optimum: %g\n", opt);
    4335
    4436    if (sense > 0) {if (opt > bound + OBBT_EPS)    {bound = (isint ? ceil (opt) : opt); return true;}}
     
    5547                OsiCuts &cs,
    5648                t_chg_bounds *chg_bds,
     49                const CoinWarmStart *warmstart,
    5750                Bonmin::BabInfo * babInfo,
    5851                int sense) {
     
    7972  // for all (original+auxiliary) variables x_i,
    8073  ////////////////////////////////////////////////////
     74
    8175  for (int i=0; i<ncols; i++)
    8276    //  for (int i=0; i<problem_ ->nVars(); i++)
     
    116110      csi -> writeLp (fname);
    117111#endif
     112
     113      csi -> setWarmStart (warmstart);
     114
    118115
    119116      if (updateBound (csi, sense, bound, isInt)) {
     
    170167
    171168      // if we solved the problem on the objective function's
    172       // auxiliary variable, it is worth updating the current point
    173       // (it will be used later to generate new cuts).
     169      // auxiliary variable (that is, we re-solved the extended
     170      // problem), it is worth updating the current point (it will be
     171      // used later to generate new cuts).
    174172      if ((objind == i) && (csi -> isProvenOptimal ()))
    175173        p -> update (const_cast <CouNumber *> (csi -> getColSolution ()), NULL, NULL);
     
    188186                               t_chg_bounds *chg_bds,
    189187                               Bonmin::BabInfo * babInfo) const {
    190   int nImprov;
    191 
     188
     189  // set large bounds to infinity (as per suggestion by JJF)
     190
     191  int ncols = csi -> getNumCols ();
     192  const double *lb = csi -> getColLower (),
     193               *ub = csi -> getColUpper ();
     194
     195  while (ncols--) {
     196    if (lb [ncols] < - COUENNE_INFINITY) csi -> setColLower (ncols, -COIN_DBL_MAX);
     197    if (ub [ncols] >   COUENNE_INFINITY) csi -> setColUpper (ncols,  COIN_DBL_MAX);
     198  }
     199
     200  //  csi -> setHintParam (OsiDoDualInResolve, false);
     201
     202  // setup cloned interface for later use
    192203  csi -> setObjSense (1); // minimization
    193204  csi -> setIntParam (OsiMaxNumIteration, MAX_OBBT_LP_ITERATION);
     
    195206  csi -> initialSolve ();
    196207
    197   nImprov =  obbt_stage (this, csi, cs, chg_bds, babInfo,  1); // minimization
    198   nImprov += obbt_stage (this, csi, cs, chg_bds, babInfo, -1); // maximization
     208  int nImprov;
     209
     210  const CoinWarmStart *warmstart = csi -> getWarmStart ();
     211
     212  // improve each bound
     213  nImprov =  obbt_stage (this, csi, cs, chg_bds, warmstart, babInfo,  1); // lower
     214  nImprov += obbt_stage (this, csi, cs, chg_bds, warmstart, babInfo, -1); // upper
    199215
    200216  return nImprov;
  • branches/Couenne/Couenne/src/convex/operators/conv-exprMul-genCuts.cpp

    r559 r593  
    104104      yi = ye -> Index ();
    105105
     106  // if expression is x*c or c*y, with c constant from the problem
     107  // definition or from the branching rules, the expression is
     108  // linear. Add one convexification equality constraint.
     109
     110  // check if either operand is constant
     111
     112  bool is0const = (xe -> Type () == CONST),
     113       is1const = (ye -> Type () == CONST);
     114
     115  CouNumber c0, c1;
     116
     117  // is one of the two constant?
     118
     119  if (is0const) c0 = xe -> Value (); 
     120  if (is1const) c1 = ye -> Value (); 
     121
     122  // compute bounds
     123
     124  xe -> getBounds (xle, xue);
     125  ye -> getBounds (yle, yue);
     126  w  -> getBounds (wle, wue);
     127
     128  CouNumber xl = (*xle) (), xu = (*xue) (),
     129            yl = (*yle) (), yu = (*yue) (),
     130            wl = (*wle) (), wu = (*wue) ();
     131
     132  delete xle; delete xue;
     133  delete yle; delete yue;
     134  delete wle; delete wue;
     135
     136  // check if either operator got constant because of the branching
     137  // rules:
     138
     139  bool i0s, i1s = i0s = false;
     140
     141  if (!is0const) {
     142
     143    if (is1const) i0s = (fabs (c1) * (xu-xl) < COUENNE_EPS);
     144    else          i0s = ((yu-yl)   * (xu-xl) < COUENNE_EPS);
     145
     146    if (i0s)
     147      c0 = 0.5 * (xl+xu);
     148  }
     149
     150  // and y
     151
     152  if (!is1const) {
     153
     154    if (is0const) i1s = (fabs (c0) * (yu-yl) < COUENNE_EPS);
     155    else          i1s = ((xu-xl)   * (yu-yl) < COUENNE_EPS);
     156
     157    if (i1s)
     158      c1 = 0.5 * (yl+yu);
     159  }
     160
     161  if (i0s) is0const = true;
     162  if (i1s) is1const = true;
     163
     164  // right now c0 and c1 only have a value if the corresponding
     165  // expression is constant
     166
     167  if (is0const || is1const) {
     168
     169    if (cg -> isFirst () ||            // if first call or
     170        ((xe -> Type () != CONST) &&   // neither term is a defined constant
     171         (ye -> Type () != CONST))) {  // (=> implied by branching rule)
     172
     173      if (is0const && is1const)
     174
     175        // w = c0*c1, which is either because the intervals got very
     176        // narrow, or because these are indeed constant (which should
     177        // have been dealt with in simplify(), but who knows...)
     178
     179        cg -> createCut (cs, c0 * c1, 0, wi, 1.);
     180
     181      else {
     182
     183        CouNumber coe;
     184        int ind;
     185
     186        if (is0const) {coe = c0; ind = yi;} // c*y
     187        else          {coe = c1; ind = xi;} // x*c
     188
     189        cg -> createCut (cs, 0., 0, wi, -1., ind, coe);
     190      }
     191    }
     192
     193    return;
     194  }
     195
    106196  bool cLX,  cRX,  cLY,  cRY,  cLW,  cRW =
    107197       cLX = cRX = cLY = cRY = cLW = true;
     
    113203  }
    114204
    115   // if expression is x*c or c*y, with c constant from the problem
    116   // definition or from the branching rules, the expression is
    117   // linear. Add one convexification equality constraint.
    118 
    119   // check if either operand is constant
    120 
    121   bool is0const = (xe -> Type () == CONST),
    122        is1const = (ye -> Type () == CONST);
    123 
    124   CouNumber c0, c1;
    125 
    126   // is one of the two constant?
    127 
    128   if (is0const) c0 = xe -> Value (); 
    129   if (is1const) c1 = ye -> Value (); 
    130 
    131   // compute bounds
    132 
    133   xe -> getBounds (xle, xue);
    134   ye -> getBounds (yle, yue);
    135   w  -> getBounds (wle, wue);
    136 
    137   CouNumber xl = (*xle) (), xu = (*xue) (),
    138             yl = (*yle) (), yu = (*yue) (),
    139             wl = (*wle) (), wu = (*wue) ();
    140 
    141   delete xle; delete xue;
    142   delete yle; delete yue;
    143   delete wle; delete wue;
    144 
    145   // check if either operator got constant because of the branching
    146   // rules:
    147 
    148   bool i0s, i1s = i0s = false;
    149 
    150   if (!is0const) {
    151 
    152     if (is1const) i0s = (fabs (c1) * (xu-xl) < COUENNE_EPS);
    153     else          i0s = ((yu-yl)   * (xu-xl) < COUENNE_EPS);
    154 
    155     if (i0s)
    156       c0 = 0.5 * (xl+xu);
    157   }
    158 
    159   // and y
    160 
    161   if (!is1const) {
    162 
    163     if (is0const) i1s = (fabs (c0) * (yu-yl) < COUENNE_EPS);
    164     else          i1s = ((xu-xl)   * (yu-yl) < COUENNE_EPS);
    165 
    166     if (i1s)
    167       c1 = 0.5 * (yl+yu);
    168   }
    169 
    170   if (i0s) is0const = true;
    171   if (i1s) is1const = true;
    172 
    173   // right now c0 and c1 only have a value if the corresponding
    174   // expression is constant
    175 
    176   if (is0const || is1const) {
    177 
    178     if (cg -> isFirst () ||            // if first call or
    179         ((xe -> Type () != CONST) &&   // neither term is a defined constant
    180          (ye -> Type () != CONST))) {  // (=> implied by branching rule)
    181 
    182       if (is0const && is1const)
    183 
    184         // w = c0*c1, which is either because the intervals got very
    185         // narrow, or because these are indeed constant (which should
    186         // have been dealt with in simplify(), but who knows...)
    187 
    188         cg -> createCut (cs, c0 * c1, 0, wi, 1.);
    189 
    190       else {
    191 
    192         CouNumber coe;
    193         int ind;
    194 
    195         if (is0const) {coe = c0; ind = yi;} // c*y
    196         else          {coe = c1; ind = xi;} // x*c
    197 
    198         cg -> createCut (cs, 0., 0, wi, -1., ind, coe);
    199       }
    200     }
    201 
    202     return;
    203   }
    204 
    205205  // Add McCormick convexification cuts:
    206206  //
     
    214214
    215215  if ((cLX || cLY) && is_boundbox_regular (yl, xl)) cg -> createCut (cs, yl*xl,-1,wi,-1.,xi,yl,yi,xl);
    216   if ((cLX || cLY) && is_boundbox_regular (yu, xu)) cg -> createCut (cs, yu*xu,-1,wi,-1.,xi,yu,yi,xu);
    217   if ((cLX || cLY) && is_boundbox_regular (yl, xu)) cg -> createCut (cs, yl*xu,+1,wi,-1.,xi,yl,yi,xu);
    218   if ((cLX || cLY) && is_boundbox_regular (yu, xl)) cg -> createCut (cs, yu*xl,+1,wi,-1.,xi,yu,yi,xl);
     216  if ((cRX || cRY) && is_boundbox_regular (yu, xu)) cg -> createCut (cs, yu*xu,-1,wi,-1.,xi,yu,yi,xu);
     217  if ((cRX || cLY) && is_boundbox_regular (yl, xu)) cg -> createCut (cs, yl*xu,+1,wi,-1.,xi,yl,yi,xu);
     218  if ((cLX || cRY) && is_boundbox_regular (yu, xl)) cg -> createCut (cs, yu*xl,+1,wi,-1.,xi,yu,yi,xl);
    219219
    220220  // add different cuts, to cut out current point in bounding box but
Note: See TracChangeset for help on using the changeset viewer.