impliedBounds-exprSum.cpp

Go to the documentation of this file.

/* $Id: impliedBounds-exprSum.cpp 942 2013-02-02 23:40:27Z pbelotti $
 *
 * Name:    impliedBounds-exprSum.cpp
 * Author:  Pietro Belotti
 * Purpose: implied bound enforcing for exprSum and exprGroup
 *
 * (C) Carnegie-Mellon University, 2006-10.
 * This file is licensed under the Eclipse Public License (EPL)
 */

#include "CouenneExprSum.hpp"
#include "CouenneExprGroup.hpp"
#include "CouenneConfig.h"
#include "CoinHelperFunctions.hpp"
#include "CoinFinite.hpp"

using namespace Couenne;

static CouNumber scanBounds (int, int, int *, CouNumber *, CouNumber *, int *);



bool exprSum::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg, enum auxSign sign) {

  // compute number of pos/neg coefficients and sum up constants

  CouNumber
    a0 = 0.,   // constant term in the sum
    wl = sign == expression::AUX_GEQ ? -COIN_DBL_MAX : l [wind],
    wu = sign == expression::AUX_LEQ ?  COIN_DBL_MAX : u [wind];

  // quick check: if both are infinite, nothing is implied...

  if ((wl < -COUENNE_INFINITY) &&
      (wu >  COUENNE_INFINITY))
    return false;

  int 
    nterms = nargs_, // # nonconstant terms
    nlin   = 0;      // # linear terms

  exprGroup *eg = NULL;

  bool isExprGroup = (code () == COU_EXPRGROUP);

  if (isExprGroup) { // if exprGroup, compute no. linear terms

    eg = dynamic_cast <exprGroup *> (this);

    a0   += eg -> getc0 ();
    nlin += eg -> lcoeff (). size ();
  }

  nterms += nlin;

  // Coefficients and indices of the positive and the negative
  // non-constant parts of the sum (at most nlin are negative, as all
  // "nonlinear" terms have coefficient 1)

  CouNumber *C1 = (CouNumber *) malloc (nterms * sizeof (CouNumber)),
            *C2 = (CouNumber *) malloc (nlin   * sizeof (CouNumber));
  int       *I1 = (int       *) malloc (nterms * sizeof (int)),
            *I2 = (int       *) malloc (nlin   * sizeof (int));

  int ipos, ineg = ipos = 0; // #pos and #neg terms

  std::set <int> intSet; // contains indices of integer variables

  // classify terms as positive or constant for the exprSum

  for (int i = nargs_; i--;) {

    int index = arglist_ [i] -> Index ();

    if (index < 0) // assume a non variable is a constant
      a0 += arglist_ [i] -> Value ();
    else {
      I1 [ipos]   = index;
      C1 [ipos++] = 1.;

      // add entry to integer variable
      if (arglist_ [i] -> isDefinedInteger ())
        intSet.insert (index);
    }
  }

  // classify terms as pos/neg or constant for the remaining of exprGroup

  if (isExprGroup) {

    exprGroup::lincoeff &lcoe = eg -> lcoeff ();

    for (register exprGroup::lincoeff::iterator el = lcoe.begin ();
         el != lcoe.end (); ++el) {

      register CouNumber coe = el -> second;
      register int       ind = el -> first -> Index ();

      if      (coe >  0.) {I1 [ipos] = ind; C1 [ipos++] = coe;}
      else if (coe < -0.) {I2 [ineg] = ind; C2 [ineg++] = coe;}

      if (el -> first -> isDefinedInteger ())
        intSet.insert (ind);
    }
  }

  // realloc to save some memory

  /*printf ("  a0 = %g, w in [%g, %g]\n", a0, wl, wu);
  printf ("  pos: "); for (int i=0; i<ipos; i++) printf ("%g x%d [%g,%g] ", C1 [i], I1 [i], l [I1 [i]], u [I1 [i]]); printf ("\n");
  printf ("  neg: "); for (int i=0; i<ineg; i++) printf ("%g x%d [%g,%g] ", C2 [i], I2 [i], l [I2 [i]], u [I2 [i]]); printf ("\n");*/

  // now we have correct  I1, I2, C1, C2, ipos, ineg, and a0

  // indices of the variable in I1 or I2 with infinite lower or upper
  // bound. If more than one is found, it is set to -2

  int infLo1 = -1, infLo2 = -1,
      infUp1 = -1, infUp2 = -1;

  // Upper bound of the sum, considering lower/upper bounds of the
  // variables but negliging the infinite ones:
  //
  // lower = $a0 + \sum_{i in I_1: a_i <  \infinity} a_i l_i
  //             + \sum_{i in I_2: a_i > -\infinity} a_i u_i$
  //
  // upper = $a0 + \sum_{i in I_1: a_i <  \infinity} a_i u_i
  //             + \sum_{i in I_2: a_i > -\infinity} a_i l_i$

  CouNumber

    lower = a0 + scanBounds (ipos, -1, I1, C1, l, &infLo1)
               + scanBounds (ineg, +1, I2, C2, u, &infUp2),

    upper = a0 + scanBounds (ipos, +1, I1, C1, u, &infUp1)
               + scanBounds (ineg, -1, I2, C2, l, &infLo2);

  // Compute lower bound for all or for some of the variables: There
  // is a bound for all variables only if both infUp1 and infLo2 are
  // -1, otherwise there is a bound only for one variable if one is -1
  // and the other is nonnegative.

  bool tighter = false;

  // check if there is room for improvement

  {
    CouNumber 
      slackL = lower - wl, // each is positive if no implication 
      slackU = wu - upper; // 

    // if lower < wl or upper > wu, some bounds can be tightened.
    // 
    // otherwise, there is no implication, but if lower
    //
    // steal some work to bound propagation... 

    if ((slackL > 0.) &&
        (infLo1 == -1) && 
        (infUp2 == -1)) {   // no implication on lower

      // propagate lower bound to w
      if ((sign != expression::AUX_LEQ) && (updateBound (-1, l + wind, lower))) {
        chg [wind].setLower(t_chg_bounds::CHANGED);
        tighter = true;
        if (intSet.find (wind)!= intSet.end ()) 
          l [wind] = ceil (l [wind] - COUENNE_EPS);
      }

      if ((slackU > 0.) &&
          (infLo2 == -1) && 
          (infUp1 == -1)) { // no implication on upper

        // propagate upper bound to w
        if ((sign != expression::AUX_GEQ) && (updateBound (+1, u + wind, upper))) {
          chg [wind].setUpper(t_chg_bounds::CHANGED);
          tighter = true;
          if (intSet.find (wind)!= intSet.end ()) 
            u [wind] = floor (u [wind] + COUENNE_EPS);
        }

        free (I1); free (I2);
        free (C1); free (C2);

        return false; // both bounds were weak, no implication possible
      }
    }
    else if ((slackU > 0.) && 
             (infLo2 == -1) && 
             (infUp1 == -1)) {

      // propagate upper bound to w
      if ((sign != expression::AUX_GEQ) && (updateBound (+1, u + wind, upper))) {
        tighter = true;
        chg [wind].setUpper(t_chg_bounds::CHANGED);
        if (intSet.find (wind)!= intSet.end ()) 
          u [wind] = floor (u [wind] + COUENNE_EPS);
      }
    }
  }

  // Subtle... make two copies of lower and upper to avoid updating
  // bounds with some previously updated bounds

  // first of all, find maximum index in I1 and I2

  int maxind = -1;

  for (register int i=ipos; i--; I1++) if (*I1 > maxind) maxind = *I1;
  for (register int i=ineg; i--; I2++) if (*I2 > maxind) maxind = *I2;

  I1 -= ipos;
  I2 -= ineg;

  CouNumber *lc = (CouNumber *) malloc (++maxind * sizeof (CouNumber));
  CouNumber *uc = (CouNumber *) malloc (maxind   * sizeof (CouNumber));

  CoinCopyN (l, maxind, lc);
  CoinCopyN (u, maxind, uc);

  // Update lowers in I1 and uppers in I2

  if ((infLo1 == -1) && (infUp2 == -1) && (wu < COUENNE_INFINITY / 1e10)) { 
    // All finite bounds. All var. bounds can be tightened.

    // tighten upper bound of variables in I1
    for (register int i=ipos; i--;) {
      int ind = I1 [i];
      if ((tighter = (updateBound (+1, u + ind, (wu - lower) / C1 [i] + lc [ind]) || tighter))) {
        chg [ind].setUpper(t_chg_bounds::CHANGED);
        if (intSet.find (ind)!= intSet.end ()) 
          u [ind] = floor (u [ind] + COUENNE_EPS);
      }
    }

    // tighten lower bound of variables in I2
    for (register int i=ineg; i--;) {
      int ind = I2 [i];
      if ((tighter = (updateBound (-1, l + ind, (wu - lower) / C2 [i] + uc [ind]) || tighter))) {
        chg [ind].setLower(t_chg_bounds::CHANGED);
        if (intSet.find (ind)!= intSet.end ()) 
          l [ind] = ceil (l [ind] - COUENNE_EPS);
      }
    }
  } else

    if ((infLo1 >= 0) && (infUp2 == -1)) {    // There is one infinite lower bound in I1
      int ind = I1 [infLo1];
      if ((tighter = (updateBound (+1, u + ind, (wu - lower) / C1 [infLo1]) || tighter))) {
        chg [ind].setUpper(t_chg_bounds::CHANGED);
        if (intSet.find (ind)!= intSet.end ()) 
          u [ind] = floor (u [ind] + COUENNE_EPS);
      }
    }
    else 
      if ((infLo1 == -1) && (infUp2 >= 0)) {  // There is one infinite upper bound in I2
        int ind = I2 [infUp2];
        if ((tighter = (updateBound (-1, l + ind, (wu - lower) / C2 [infUp2]) || tighter))) {
          chg [ind].setLower(t_chg_bounds::CHANGED);
          if (intSet.find (ind)!= intSet.end ()) 
            l [ind] = ceil (l [ind] - COUENNE_EPS);
        }
      }

  // Update uppers in I1 and lowers in I2

  if ((infUp1 == -1) && (infLo2 == -1) && (wl > -COUENNE_INFINITY / 1e10)) { 
    // All finite bounds. All var. bounds can be tightened.

    for (register int i=ipos; i--;) {
      int ind = I1 [i];
      if ((tighter = (updateBound (-1, l + ind, (wl - upper) / C1 [i] + uc [ind]) || tighter))) {
        chg [ind].setLower(t_chg_bounds::CHANGED); 
        if (intSet.find (ind) != intSet.end ()) 
          l [ind] = ceil (l [ind] - COUENNE_EPS);
      }
    }

    for (register int i=ineg; i--;) {
      int ind = I2 [i];
      if ((tighter = (updateBound (+1, u + ind, (wl - upper) / C2 [i] + lc [ind]) || tighter))) {
        chg [ind].setUpper(t_chg_bounds::CHANGED);
        if (intSet.find (ind) != intSet.end ()) 
          u [ind] = floor (u [ind] + COUENNE_EPS);
      }
    }
  } else 

    if ((infUp1 >= 0) && (infLo2 == -1)) { // There is one infinite lower bound in I2
      int ind = I1 [infUp1];
      if ((tighter = (updateBound (-1, l + ind, (wl - upper) / C1 [infUp1]) || tighter))) {
        chg [ind].setLower(t_chg_bounds::CHANGED);
        if (intSet.find (ind) != intSet.end ()) 
          l [ind] = ceil (l [ind] - COUENNE_EPS);
      }
    }
    else 
      if ((infUp1 == -1) && (infLo2 >= 0)) {  // There is one infinite upper bound in I1
        int ind = I2 [infLo2];
        if ((tighter = (updateBound (+1, u + ind, (wl - upper) / C2 [infLo2]) || tighter))) {
          chg [ind].setUpper(t_chg_bounds::CHANGED);
          if (intSet.find (ind) != intSet.end ()) 
            u [ind] = floor (u [ind] + COUENNE_EPS);
        }
      }

  // ...phew!

  free (I1); free (I2);
  free (C1); free (C2);
  free (lc); free (uc);

  return tighter;
}



static CouNumber scanBounds (int        num,      
                             int        sign,     
                             int       *indices,  
                             CouNumber *coeff,    
                             CouNumber *bounds,   
                             int       *infnum) { 
  CouNumber bound = 0.;

  for (register int i = num; i--;) {

    CouNumber bd = bounds [indices [i]];

    // be sensitive here, check for bounds a little within the finite realm

    if (((sign > 0) ? bd : -bd) > COUENNE_INFINITY / 1e10) {

      bounds [indices [i]] = (sign > 0) ? COIN_DBL_MAX : -COIN_DBL_MAX;
      //bounds [indices [i]] = (sign > 0) ? 1e300 : -1e300;

      // this variable has an infinite bound, mark it
      if      (*infnum == -1) *infnum =  i; // first variable with infinite bound, so far
      else if (*infnum >=  0) *infnum = -2; // oops... more than one found, no finite bound
    }
    else bound += coeff [i] * bd; // no infinity detected, sum a weighted, finite bound
  }

  return bound;
}