exprPow.cpp

Go to the documentation of this file.

/* $Id: exprPow.cpp 1147 2015-05-04 14:01:51Z stefan $
 *
 * Name:    exprPow.cpp
 * Author:  Pietro Belotti
 * Purpose: definition of powers
 *
 * (C) Carnegie-Mellon University, 2006-11.
 * This file is licensed under the Eclipse Public License (EPL)
 */

#include <math.h>
#include <assert.h>

#include "CouennePrecisions.hpp"
#include "CouenneExprPow.hpp"
#include "CouenneExprSum.hpp"
#include "CouenneExprMul.hpp"
#include "CouenneExprDiv.hpp"
#include "CouenneExprLog.hpp"
#include "CouenneExprConst.hpp"
#include "CouenneProblem.hpp"

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

using namespace Couenne;


expression *exprPow::simplify () {

  exprOp:: simplify ();

  if ((*arglist_) -> Type () == CONST) { // expr = c0 ^ g(x)

    CouNumber c0 = (*arglist_) -> Value ();

    if (arglist_ [1] -> Type () == CONST) { // expr = c0 ^ c1

      CouNumber c1 = arglist_ [1] -> Value ();

      // delete arglist_ [0]; 
      // delete arglist_ [1];

      // arglist_ [0] = arglist_ [1] = NULL;

      return new exprConst
        (issignpower_ ? 
         COUENNE_sign(c0) * pow (fabs(c0), c1) : 
         pow (c0, c1));
    }
    else 
      if (fabs (c0) <= COUENNE_EPS_SIMPL) 
        return new exprConst (0.);
  }
  else // only need to check if g(x) == 0

    if (arglist_ [1] -> Type () == CONST) {

      CouNumber expon = arglist_ [1] -> Value ();

      if (fabs (expon) <= COUENNE_EPS_SIMPL) // expr = x ^ 0 = 1
        return new exprConst (1.);

      else if (fabs (expon - 1.) <= COUENNE_EPS_SIMPL) { // expr = x ^ 1 = x

        //delete arglist_ [1];
        expression *ret = arglist_ [0];
        arglist_ [0] = NULL;

        return ret;
      }

      else if (fabs (expon + 1.) <= COUENNE_EPS_SIMPL) { // expr = x ^ -1 = 1/x

        //delete arglist_ [1];
        expression *ret = new exprInv (arglist_ [0]);
        arglist_ [0] = NULL;
        return ret;
      }

      //
      // x^k = x for x binary. Too bad we don't know bounds yet, so the code below will give segfault

      // TODO: readnl first reads bounds. Then uncomment this.

      //       // is it an integer variable with bounds [-1,0] or [0,1]
      //       else if ((arglist_ [0] -> Type () == VAR) && (arglist_ [0] -> isDefinedInteger ())) {

      //        CouNumber lb, ub;
      //        arglist_ [0] -> getBounds (lb, ub);

      //        if ((fabs (lb)      < COUENNE_EPS) &&
      //            (fabs (ub - 1.) < COUENNE_EPS)) {  // {0,1}

      //          delete arglist_ [1];
      //          expression *ret = arglist_ [0];
      //          arglist_ [0] = arglist_ [1] = NULL;
      //          return ret;

      //        } else if ((fabs (lb + 1.) < COUENNE_EPS) &&
      //                   (fabs (ub)      < COUENNE_EPS)) { // {-1,0}

      //          delete arglist_ [1];
      //          expression *ret = new exprOpp (arglist_ [0]);
      //          arglist_ [0] = arglist_ [1] = NULL;
      //          return ret;
      //        }
      //       }

    }

  return NULL;
}



expression *exprPow::differentiate (int index) {

  if (!(arglist_ [0] -> dependsOn (index))  &&
      !(arglist_ [1] -> dependsOn (index)))
    return new exprConst (0.);

  if (arglist_ [0] -> Type () == CONST) { // k^f(x), k constant

    CouNumber base = arglist_ [0] -> Value ();

    if (base == 0.)
      return new exprConst (0.);

    return new exprMul (new exprConst (log (base)),
                        new exprMul (new exprPow (new exprConst (base), 
                                                  arglist_ [1] -> clone ()),
                                     arglist_ [1] -> differentiate (index)));

  } else if (arglist_ [1] -> Type () == CONST) { // f(x)^k, k constant

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

    return new exprMul (new exprConst (exponent),
                        new exprMul (new exprPow (arglist_ [0] -> clone (),
                                                  new exprConst (exponent - 1.)),
                                     arglist_ [0] -> differentiate (index)));
  }

  // all other cases: f(x)^g(x)

  expression **alm  = new expression * [2];
  expression **alp  = new expression * [2];
  expression **als  = new expression * [2];
  expression **alm1 = new expression * [2];
  expression **alm2 = new expression * [2];
  expression **ald  = new expression * [2];

  alp [0] = new exprClone (arglist_ [0]);
  alp [1] = new exprClone (arglist_ [1]);

  alm [0] = new exprPow (alp, 2);

  alm1 [0] = arglist_ [1] -> differentiate (index);
  alm1 [1] = new exprLog (new exprClone (arglist_ [0]));

  als [0] = new exprMul (alm1, 2);

  ald [0] = new exprClone (arglist_ [1]);
  ald [1] = new exprClone (arglist_ [0]);

  alm2 [0] = new exprDiv (ald, 2);
  alm2 [1] = arglist_ [0] -> differentiate (index);

  als [1] = new exprMul (alm2, 2);

  alm [1] = new exprSum (als, 2);

  return new exprMul (alm, 2);
}



int exprPow::Linearity () {

  if (arglist_ [0] -> Type () == CONST) {

    if (arglist_ [1] -> Type () == CONST) return CONSTANT;
    else                                  return NONLINEAR;
  }
  else {

    double exponent = arglist_ [1] -> Value ();

    if (fabs (exponent - COUENNE_round (exponent)) > COUENNE_EPS)
      return NONLINEAR;

    if (arglist_ [1] -> Type () == CONST) { 

      int expInt = (int) COUENNE_round (exponent);

      if (arglist_ [0] -> Linearity () == LINEAR) {

        switch (expInt) {

        case 0:  return CONSTANT;
        case 1:  return LINEAR;
        case 2:  return (issignpower_ ? NONLINEAR : QUADRATIC);

        default: return NONLINEAR;
        }
      }
      else 
        if (arglist_ [0] -> Linearity () == QUADRATIC) 
          switch (expInt) {

          case 0:  return CONSTANT;
          case 1:  return QUADRATIC;

          default: return NONLINEAR;
          }
        else return NONLINEAR;
    }
    else return NONLINEAR;
  }
}


bool exprPow::isInteger () {

  // base

  if (!(arglist_ [0] -> isInteger ())) { 

    // base not integer: check if constant and integer
    CouNumber lb, ub;
    arglist_ [0] -> getBounds (lb, ub);

    if ((fabs (lb - ub) > COUENNE_EPS) ||
        !::isInteger (lb))
      return false;
  }

  // exponent

  if (!(arglist_ [1] -> isInteger ())) { 

    // exponent not defined integer: check if constant and at integer
    // value (and positive, or base negative integer)

    CouNumber lb, ub;
    arglist_ [1] -> getBounds (lb, ub);

    if ((fabs (lb - ub) > COUENNE_EPS) ||
        !::isInteger (lb))
      return false;

    if (lb < 0) { // exponent negative, check again base

      arglist_ [0] -> getBounds (lb, ub);

      if ((fabs (lb - ub) > COUENNE_EPS) ||
          (fabs (lb) < COUENNE_EPS) ||
          !::isInteger (1. / lb))
        return false;
    }
  } else {

    // if base integer and exponent integer, must check that exponent
    // is nonnegative

    CouNumber lb, ub;
    arglist_ [1] -> getBounds (lb, ub);

    if (lb < .5)
      return false;
  }

  return true;
}


void exprPow::closestFeasible (expression *varind,
                               expression *vardep, 
                               CouNumber &left,
                               CouNumber &right) const {
  CouNumber
    x  = (*varind) (),
    y  = (*vardep) (),
    k  = arglist_ [1] -> Value (),
    xk = safe_pow (x,    k, issignpower_),
    yk = safe_pow (y, 1./k, issignpower_);

  assert(!issignpower_ || k > 0);

  int intk = 0;

  bool isInt    =            fabs (k    - (double) (intk = COUENNE_round (k)))    < COUENNE_EPS,
       isInvInt = !isInt && (fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS);

  // three cases: 
  // 1) k or  1/k odd or signpower => have either left or right
  // 2) k or  1/k even,            => may have both
  // 3) k and 1/k fractional       => have either left or right

  if (isInt || isInvInt)

    if (intk % 2 || issignpower_) // case 1

      if (k > 0) 
        ((y < xk) ? left : right) = yk; // easy, x^k is continuous

      else

        if      (y < 0.)          // third, fourth orthant
          if (y < xk) right = yk; // in convex region y < 1/x within third orthant
          else        left  = yk; // remaining non-convex area

        else                      // first, second orthant
          if (y > xk) left  = yk; // in convex region y > 1/x within first orthant
          else        right = yk; // remaining non-convex area

    else // case 2

      if (y <= 0.) // third, fourth orthant => no solution
        left = - (right = COIN_DBL_MAX);

      else

        if (k > 0) 

          if (k < 1) // roots, have x >= 0

            if (x > yk) left  = yk;
            else        right = yk;

          else

            if (x > yk)       left  =  yk;
            else if (x < -yk) right = -yk;
            else              left  = - (right = yk);

        else // k negative
          if (y < xk) // between asymptotes
            left = - (right = yk);
          else  // in one of the two convex areas
            if (x > 0) left  =  yk;
            else       right = -yk;

  else // case 3: assume x bounded from below by 0

    if (k > 0) ((y < xk) ? left : right) = yk;
    else       ((y > xk) ? left : right) = yk;
}


CouNumber exprPow::gradientNorm (const double *x) {

  int ind0 = arglist_ [0] -> Index ();
  CouNumber exponent = arglist_ [1] -> Value ();
  return (ind0 < 0) ? 0. : fabs (exponent * safe_pow (x [ind0], exponent - 1, issignpower_));
}


bool exprPow::isCuttable (CouenneProblem *problem, int index) const {

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

  bool
    isInt    = ::isInteger (exponent),
    isInvInt = (exponent != 0.) && ::isInteger (1. / exponent);

  int intExp = (isInt ? COUENNE_round (exponent) : (isInvInt ? COUENNE_round (1. / exponent) : 0));

  if (exponent > 0.) {

    if (isInt || isInvInt) {

      // TODO also with an odd exponent, there are convex parts where one could do separation
      if (intExp % 2 || issignpower_) return false; // exponent odd or 1/odd

      CouNumber 
        x = problem -> X (arglist_ [0] -> Index ()),
        y = problem -> X (index);

      if (isInt) return (y <= safe_pow (x, exponent, issignpower_)); // below convex curve ==> cuttable

      return (y >= safe_pow (x, exponent, issignpower_)); // above concave k-th root curve ==> cuttable
    } else {

      // non-integer exponent
      CouNumber 
        x = problem -> X (arglist_ [0] -> Index ()),
        y = problem -> X (index);

      return (((exponent <= 1.) && (y >= safe_pow (x, exponent))) ||
              ((exponent >= 1.) && (y <= safe_pow (x, exponent))));
    }      
  } else {

    // non-integer exponent
    CouNumber 
      x  = problem -> X (arglist_ [0] -> Index ()),
      y  = problem -> X (index),
      lb = problem -> Lb (index),
      ub = problem -> Ub (index);

    if (isInt || isInvInt)

      if (!(intExp % 2 || issignpower_)) return (((lb > 0) || (ub < 0)) && (y * safe_pow (fabs (x), -exponent, issignpower_) <= 1.));
      else               return (((lb > 0) || (ub < 0)) && (y * safe_pow (x,        -exponent, issignpower_) <= 1.));
    else                 return                            (y * safe_pow (x,        -exponent, issignpower_) <= 1.);
  }
}