impliedBounds-exprPow.cpp

Go to the documentation of this file.

/* $Id: impliedBounds-exprPow.cpp 1047 2014-01-20 00:22:31Z pbelotti $
 *
 * Name:    impliedBounds-exprPow.cpp
 * Author:  Pietro Belotti
 * Purpose: implied bounds for power operators
 *
 * (C) Carnegie-Mellon University, 2008-10.
 * This file is licensed under the Eclipse Public License (EPL)
 */

#include "CouenneExprPow.hpp"
#include "CouenneExpression.hpp"
#include "CouenneConfig.h"
#include "CoinFinite.hpp"
#include "CoinHelperFunctions.hpp"

using namespace Couenne;


void invPowImplBounds (int, int, CouNumber *, CouNumber *, CouNumber, bool &, bool &, enum expression::auxSign);



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

  //int xi = arglist_ [0] -> Index ();
  //if (xi>=0) printf ("in implBound-pow: %g,%g\n", l [xi], u [xi]);

  bool resL, resU = resL = false;

  if (arglist_ [0] -> Type () == CONST)   // base is constant or zero, nothing to do
    return false;

  assert (arglist_ [1] -> Type () == CONST);

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

  CouNumber k = arglist_ [1] -> Value (); // exponent

  if ((fabs (k) < COUENNE_EPS) || 
      (fabs (k) > COUENNE_INFINITY)) // a null or infinite k is of little use
    return false;

  int intk; // integer (or integer inverse of) exponent

  bool 
    isint    =           (            fabs (k    - (intk = COUENNE_round (k)))    < COUENNE_EPS),  // k   integer
    isinvint = !isint && (k != 0. && (fabs (1./k - (intk = COUENNE_round (1./k))) < COUENNE_EPS)); // 1/k integer

  CouNumber 
    wl = ((sign == expression::AUX_GEQ) ? -COIN_DBL_MAX : l [wind]), // lower w
    wu = ((sign == expression::AUX_LEQ) ?  COIN_DBL_MAX : u [wind]); // upper w

  if ((isint || isinvint) && (intk % 2 || issignpower_)) { 

    // k or 1/k integer and odd, or non-integer --> it is a monotone
    // increasing function, apart when k negative

    if (k > 0.) { // simple, just follow bounds

      if (wl > - COUENNE_INFINITY) resL = updateBound (-1, l + index, safe_pow (wl, 1./k, issignpower_)); 
      if (wu <   COUENNE_INFINITY) resU = updateBound (+1, u + index, safe_pow (wu, 1./k, issignpower_));

    } else {// slightly more complicated, resort to same method as in exprInv
      invPowImplBounds (wind, index, l, u, 1./k, resL, resU, sign);
      assert (!issignpower_);
    }
  } 
  else 
    if (isint) { // x^k, k integer and even --> non monotone

      CouNumber bound = (k<0) ? wl : wu;

      // |x| <= b^(1/k), where b is wl or wu depending on k negative
      // or positive, respectively

      if (bound > COUENNE_EPS) {

        if (fabs (bound) < COUENNE_INFINITY) {
          resL = updateBound (-1, l + index, - safe_pow (bound, 1./k, issignpower_));
          resU = updateBound (+1, u + index,   safe_pow (bound, 1./k, issignpower_));
        } /*else {
          resL = updateBound (-1, l + index, - COUENNE_INFINITY);
          resU = updateBound (+1, u + index,   COUENNE_INFINITY);
          }*/
      }

      // invert check, if bounds on x do not contain 0 we may improve them

      bound = (k>0) ? wl : wu;

      if (bound > 0.) {

        CouNumber
          xl = l [index], 
          xu = u [index],
          xb = safe_pow (bound, 1./k, issignpower_);

        if      (xl > - xb + COUENNE_EPS) resL = updateBound (-1, l + index,   xb) || resL;
        else if (xu <   xb - COUENNE_EPS) resU = updateBound ( 1, u + index, - xb) || resU;
      }

    } else { 

      // Two cases:
      // 1) x^k, k=(1/h), h integer and even, or 
      // 2) x^k, neither k nor 1/k integer

      CouNumber lb = wl, ub = wu;

      if (k < 0) { // swap bounds as they swap on the curve x^k when 
        lb = wu;
        ub = wl;
      }

      if ((lb > 0. || k > 0.) && (sign != expression::AUX_GEQ)) 
        resL = updateBound (-1, l + index, safe_pow (lb, 1./k, issignpower_));

      if ((fabs (ub) < COUENNE_INFINITY) && 
          (ub > 0 || k > 0.) &&
          (sign != expression::AUX_LEQ)) 
        resU = updateBound (+1, u + index, safe_pow (ub, 1./k, issignpower_));
      //else                  resU = updateBound (+1, u + index, COUENNE_INFINITY);
    }

  if (resL) chg [index].setLower(t_chg_bounds::CHANGED);
  if (resU) chg [index].setUpper(t_chg_bounds::CHANGED);

  bool xInt = arglist_ [0] -> isInteger ();

  if ((resL || resU) && xInt) {

    // careful with what "integer" means when a bound is 1e-8 (see minlplib/deb[789].nl)
    if (resL && (fabs (l [index]) > COUENNE_EPS)) l [index] = ceil  (l [index] - COUENNE_EPS);
    if (resU && (fabs (u [index]) > COUENNE_EPS)) u [index] = floor (u [index] + COUENNE_EPS);
  }

  return (resL || resU);
}