TwoImpliedIneqCombine.cpp

Go to the documentation of this file.

/* $Id: TwoImpliedIneqCombine.cpp 939 2013-01-12 19:53:30Z pbelotti $
 *
 * Name:    TwoImpliedCombine.cpp
 * Author:  Pietro Belotti
 * Purpose: Bound reduction using two inequalities from the LP relaxation
 * 
 * (C) Pietro Belotti, 2010.
 * This file is licensed under the Eclipse Public License (EPL)
 */

#include <stdlib.h>

#include "CglCutGenerator.hpp"
#include "CouenneTwoImplied.hpp"
#include "CoinPackedMatrix.hpp"
#include "CouennePrecisions.hpp"
#include "CouenneProblem.hpp"
#include "CouenneExprVar.hpp"

#define MIN_DENOM 1.e-10

namespace Couenne {

enum signum {NEG, POS, DN, UP};

//#define DEBUG

// type for alphas ////////////////////////////////////////////////

typedef struct {

  double alpha;
  int indVar;
  enum signum sign; // defines behavior of convex combination within
                    // [0,1]: always negative, always positive, + ->
                    // -, and - -> +
} threshold;


// compare two threshold values ///////////////////////////////////

int compthres (register const void *t1, 
               register const void *t2) {

  register double
    a1 = (*(threshold **) t1) -> alpha,
    a2 = (*(threshold **) t2) -> alpha;

  return ((a1 < a2) ? -1 : 
          (a1 > a2) ?  1 : 0);
}

// structure to sort indices with respective values

struct indPosPair {

  int index;
  int position;
};

// compare two pairs (index,value) ///////////////////////////////////

int compPair (register const void *p1, 
              register const void *p2) {

  register int
    i1 = ((struct indPosPair *) p1) -> index,
    i2 = ((struct indPosPair *) p2) -> index;

  return ((i1 < i2) ? -1 : 
          (i1 > i2) ?  1 : 0);
}


//
// Combines two inequalities whose indices, bounds, and coefficients
// are defined by the arguments ind, l and u, and a, respectively.
//

int combine (CouenneProblem *p,
             int n1, 
             int n2, 
             const int *ind1c, // indices
             const int *ind2c, 
             double *sa1, // coeff (sparse array) 
             double *sa2, // DO NOT INVERT (already done at caller)
             const double *a1c,  // coeff
             const double *a2c, 
             double *clb, // variable bounds
             double *cub,
             double l1, // constraint bounds
             double l2,  
             double u1, 
             double u2,
             bool *isInteger,
             int sign) { // invert second constraint? -1: yes, +1: no

  // first, sort ind1/a1 and ind2/a2 w.r.t. indices. They may not be
  // sorted and this messes up the while loop below.

  int 
    *ind1 = new int [n1],
    *ind2 = new int [n2];

  double 
    *a1 = new double [n1],
    *a2 = new double [n2];

  CouNumber 
    *Lb = p -> Lb (),
    *Ub = p -> Ub ();

  struct indPosPair *pairs = new struct indPosPair [CoinMax (n1,n2)];

  // re-order ind1 and a1, ind2 and a2 ///////////////////////////////////////////

  for (int i=0; i<n1; i++) {
    pairs [i]. index    = ind1c [i];
    pairs [i]. position = i;
  }

  // This is not very wise: qsort behaves horribly (that is, O(n^2))
  // on already sorted arrays (if the pivot element is the last one).

  qsort (pairs, n1, sizeof (indPosPair), compPair);

  for (int i=0; i<n1; i++) {

    int rightpos = pairs [i]. position;

    ind1 [i] = ind1c [rightpos];
    a1   [i] = a1c   [rightpos];
  }


  for (int i=0; i<n2; i++) {
    pairs [i]. index    = ind2c [i];
    pairs [i]. position = i;
  }

  qsort (pairs, n2, sizeof (indPosPair), compPair);

  for (int i=0; i<n2; i++) {

    int rightpos = pairs [i]. position;

    ind2 [i] = ind2c [rightpos];
    a2   [i] = a2c   [rightpos];
  }

  delete [] pairs;

  // Set multiplier of second constraint, to be used with a2 but not
  // with sa2.

  double signA;

  if (sign < 0) {

    double tmp = u2;
    u2 = -l2;
    l2 = -tmp;

    signA = -1.;

  } else signA = 1.;

#ifdef DEBUG
  printf ("here are the two ineqs (sign=%d):\n", sign);
  printf (" 1: %g <=", l1);             for (int i=0; i<n1; i++) printf (" %+g x%d",         a1 [i], ind1 [i]);
  printf ("<= %g\n 2: %g <=", u1, l2);  for (int i=0; i<n2; i++) printf (" %+g x%d", signA * a2 [i], ind2 [i]); printf ("<= %g\n", u2);
#endif

  threshold
     *alphas   = new threshold  [n1 + n2], // contains all alphas (there might be at most n1+n2-1, but let's be flexible)
    **inalphas = new threshold* [n1 + n2], // points to those in [0,1[ (will be sorted)
     *curalpha = alphas;

  int 
    i1    = 0, 
    i2    = 0, 
    cnt   = 0,
    incnt = 0;

  while (i1 < n1 ||
         i2 < n2) {

    if        (i1 < n1 && (i2 == n2 || ind1 [i1] < ind2 [i2])) {

      curalpha -> alpha = 1.;
      curalpha -> sign  = (a1 [i1] < 0.) ? NEG : POS;
      curalpha -> indVar = ind1 [i1];

      cnt++;
      i1++;

    } else if (i2 < n2 && (i1 == n1 || ind2 [i2] < ind1 [i1])) {

      curalpha -> alpha = 0.;
      curalpha -> sign  = (signA * a2 [i2] < 0.) ? NEG : POS;
      curalpha -> indVar = ind2 [i2];

      cnt++;
      i2++;

    } else {

      // the two indices are equal, < n1, n2 respectively, we may have
      // an alpha in ]0,1[

      double a2i2 = signA * a2 [i2];

      curalpha -> alpha = 
        (a1 [i1] == a2i2) ? 
        -1. : 
        - a2i2 / (a1 [i1] - a2i2);

      if (curalpha -> alpha <= 0. || 
          curalpha -> alpha >= 1.)

        curalpha -> sign = (a1 [i1] > 0.) ? POS : NEG;

      else 

        curalpha -> sign = 
          (a1 [i1] == a2i2) ? 
          ((a1 [i1] < 0.) ? NEG : POS) :
          ((a1 [i1] < 0.) ? DN  : UP);

      curalpha -> indVar = ind1 [i1];

      if (curalpha -> alpha > 0. &&
          curalpha -> alpha < 1.)
        inalphas [incnt++] = curalpha;

      cnt++;
      i1++;
      i2++;
    }

#ifdef DEBUG
    printf ("(%d,%g,%s) ", 
            curalpha -> indVar, 
            curalpha -> alpha, 
            (curalpha -> sign == NEG) ? "NEG" :
            (curalpha -> sign == POS) ? "POS" :
            (curalpha -> sign == DN)  ? "DN"  : "UP");
#endif

    curalpha++;
  }

#ifdef DEBUG
  printf ("\n");
#endif

  // If none of them has an alpha in ]0,1[, nothing needs to be done.

  if (!incnt) {

    delete [] inalphas;
    delete [] alphas;
    delete [] ind1;
    delete [] ind2;
    delete [] a1;
    delete [] a2;

    return 0;
  }

  //----------------------------------------------------------------------
  // Done setting up zeros of (alpha a'_i + (1-alpha) a''_i) x_i
  //----------------------------------------------------------------------

  // Now all thresholds are defined. For all threshold in ]0,1[,
  // determine if they define a convex combination of the two
  // inequalities that allows to tighten the bounds of all variables
  // with at least one nonzero coefficient in the two ineqs.

  int ntightened = 0;

  if (incnt > 1)
    qsort (inalphas, incnt, sizeof (threshold *), compthres);

#ifdef DEBUG
  printf ("--sorted: ");
  for (int i = 0; i < incnt; i++)
    printf ("(%d,%g,%s) ", 
            inalphas [i] -> indVar, 
            inalphas [i] -> alpha, 
            (inalphas [i] -> sign == NEG) ? "NEG" :
            (inalphas [i] -> sign == POS) ? "POS" :
            (inalphas [i] -> sign == DN)  ? "DN"  : "UP");    
  printf ("\n");
#endif

  // Compute the max/min of the two constraint bodies. Consider bounds
  // according to the sign of the CONVEX COMBINATION OF THE
  // COEFFICIENT:
  //
  // a_i = alpha a1_i + (1-alpha) a2_i
  //
  // L1_min = sum {i in N} a1_i f-(i)
  // L1_max = sum {i in N} a1_i f+(i)
  //
  // L2_min = sum {i in N} a2_i f-(i)
  // L2_max = sum {i in N} a2_i f+(i)
  //
  // with f-(i) = l_i if a_i > 0 and a*_i > 0
  //                  or a_i < 0 and a*_i < 0
  //              u_i otherwise;
  //
  //      f+(i) = u_i if a_i > 0 and a*_i > 0
  //                  or a_i < 0 and a*_i < 0
  //              l_i otherwise
  //
  // and a* is a1 or a2 depending on the equation.

  double 
    minSum1 = 0., minSum2 = 0.,
    maxSum1 = 0., maxSum2 = 0.;

  int 
    mInfs = 0, // number of -inf summed to minSum 1 and 2 (both)
    pInfs = 0; //           +inf           maxSum

  // The maximum sum in the constraint's body depends on the
  // coefficients arising from the convex combination, but alpha=0
  // at initialization.

  // compute for a_1 ////////////////////////////////////////////////////////

  for (int i=0; i<n1; i++) {

    int indVar1 = ind1 [i];

    double 
      a1i  = a1  [i],
      a2i  = sa2 [indVar1],
      clb1 = clb [indVar1],
      cub1 = cub [indVar1];

    // if no corresponding term on other side, use coefficient of this
    // inequality

    if (a2i == 0.)
      a2i = a1i;

#ifdef DEBUG
    else if (fabs (a2i) < 1e-10)
      printf ("\nNumerics in a2i\n\n");
#endif

    if        (a2i < 0.) {

      bool zeroA = (fabs (sa2 [indVar1]) == 0.);

      if (cub1 >   COUENNE_INFINITY/10) {if (zeroA) mInfs ++;} else minSum1 += a1i * cub1;
      if (clb1 < - COUENNE_INFINITY/10) {if (zeroA) pInfs ++;} else maxSum1 += a1i * clb1;
                                                                                         
    } else if (a2i > 0.) {                                                               
                                                                                         
      bool zeroA = (fabs (sa2 [indVar1]) == 0.);

      if (clb1 < - COUENNE_INFINITY/10) {if (zeroA) mInfs ++;} else minSum1 += a1i * clb1;
      if (cub1 >   COUENNE_INFINITY/10) {if (zeroA) pInfs ++;} else maxSum1 += a1i * cub1;
    }
  }

  // compute for a_2 ////////////////////////////////////////////////////////

  for (int i=0; i<n2; i++) {

    int indVar2 = ind2 [i];

    double 
      a2i  = a2  [i] * signA,
      clb2 = clb [indVar2],
      cub2 = cub [indVar2];

    if (a2i < 0.) {

      if (cub2 >   COUENNE_INFINITY/10) mInfs ++; else minSum2 += a2i * cub2;
      if (clb2 < - COUENNE_INFINITY/10) pInfs ++; else maxSum2 += a2i * clb2;

    } else {

      if (clb2 < - COUENNE_INFINITY/10) mInfs ++; else minSum2 += a2i * clb2;
      if (cub2 >   COUENNE_INFINITY/10) pInfs ++; else maxSum2 += a2i * cub2;
    }
  }

  // At this point, m[in,ax]Sum2 contain the finite max/min sums of
  // the second term ax, and pInfs and mInfs the number of bounds
  // that make the min/max sum infinite. We are at alpha=0,
  // therefore this is the max/min sum of the whole convex
  // combination of the two constraint bodies.

#ifdef DEBUG
  printf ("  at alpha=zero, m[in,ax]sum[12] = (1:(%g,%g), 2:(%g,%g)), %d mInf, %d pInf\n",
          minSum1, maxSum1, 
          minSum2, maxSum2,
          mInfs, pInfs);
#endif

  // scan all alphas in ]0,1[ ///////////////////////////////////////////////

  // for all variables x[i], look for tighter lower/upper bounds on x[i]

  std::vector <std::pair <int, double> > 
    newLB, // pairs (indVar, value) of new bounds
    newUB;

  delete [] a1;
  delete [] a2;

#ifdef DEBUG
  printf ("  ");
  for (int j=0; j < n1; j++) printf ("x%d [%g,%g] ", ind1 [j], clb [ind1 [j]], cub [ind1 [j]]); printf ("--- ");
  for (int j=0; j < n2; j++) printf ("x%d [%g,%g] ", ind2 [j], clb [ind2 [j]], cub [ind2 [j]]); printf ("\n");
#endif

  for (int i = 0; i < incnt; i++) {

    threshold *inalpha = inalphas [i];

#ifdef DEBUG
    printf ("  looking at %d: (%d,%g,%s)\n", i,
            inalphas [i] -> indVar, 
            inalphas [i] -> alpha, 
            (inalphas [i] -> sign == NEG) ? "NEG" :
            (inalphas [i] -> sign == POS) ? "POS" :
            (inalphas [i] -> sign == DN)  ? "DN"  : "UP"); 

    fflush (stdout);
#endif

    // look at each inalphas, those are the only breakpoints where
    // bounds can improve. For alpha in {0,1} there is already a
    // procedure (single constraint implied bound).

    // check for improved bounds at all variables

    double alpha = inalpha -> alpha;

    // enter next interval, update dynamic max/min sums for this
    // value of alpha

    while (true) {

      inalpha = inalphas [i];

      int
        signalpha = inalpha -> sign,   // either UP or DN (not POS or NEG)
        indVar    = inalpha -> indVar; // the index of the variable whose a_i(alpha) changes sign

      double
        clbi = clb [indVar],
        cubi = cub [indVar];

      // as a variable is changing sign of (alpha a'_i + (1-alpha)
      // a''_i), so must m[in|ax]Sum[12]

      if (fabs (clbi) > 0.) {

        if (clbi < - COUENNE_INFINITY/10) {

          if (signalpha == DN) {pInfs++; mInfs--;} 
          else                 {pInfs--; mInfs++;} // sign is UP

        } else {

          double 
            sa1ic = sa1 [indVar] * clbi,
            sa2ic = sa2 [indVar] * clbi;

          if (signalpha == DN) { // means sa1 [indVar] < 0 while sa2 [indVar] > 0

            minSum1 -= sa1ic;    minSum2 -= sa2ic;
            maxSum1 += sa1ic;    maxSum2 += sa2ic;

          } else {               // sign must be UP

            minSum1 += sa1ic;    minSum2 += sa2ic;
            maxSum1 -= sa1ic;    maxSum2 -= sa2ic;
          }
        }
      }

      // same for upper bound

      if (fabs (cubi) > 0.) {

        if (cubi >   COUENNE_INFINITY/10) {

          if (signalpha == DN) {pInfs--; mInfs++;} 
          else                 {pInfs++; mInfs--;} // sign is UP

        } else {

          double 
            sa1ic = sa1 [indVar] * cubi,
            sa2ic = sa2 [indVar] * cubi;

          if (signalpha == DN) { // means sa1 [indVar] < 0 while sa2 [indVar] > 0

            minSum1 += sa1ic;    minSum2 += sa2ic;
            maxSum1 -= sa1ic;    maxSum2 -= sa2ic;

          } else { // sign must be UP

            minSum1 -= sa1ic;    minSum2 -= sa2ic;
            maxSum1 += sa1ic;    maxSum2 += sa2ic;
          }
        }
      }

      if ((i < incnt-1) && 
          (inalphas [i+1] -> alpha == alpha)) ++i;
      else break;
    }

#ifdef DEBUG
    printf ("  at alpha=%g, m[in,ax]sum[12] = (1:(%g,%g), 2:(%g,%g)), %d mInf, %d pInf\n",
            inalphas [i] -> alpha,
            minSum1, maxSum1, 
            minSum2, maxSum2,
            mInfs, pInfs);
#endif

    for (int j=0; j<n1+n2; j++) {

      // Loop on all variables that might improve a bound as a
      // result of the convex combination with the current value of
      // alpha.

      // scan all nonzeros, those on one side and those on both
      // sides. This has already been scanned.
      if (j >= n1 && sa1 [ind2 [j - n1]] != 0.) 
        continue; 

      int
        indVar = ((j < n1) ? ind1 [j] : ind2 [j - n1]),
        tickMin = 0, // one if subtracting this term removes one infinity from minSum
        tickMax = 0; //                                                        max   

      double
        sa1i = sa1 [indVar],
        sa2i = sa2 [indVar],
        ci = alpha * sa1i + (1. - alpha) * sa2i; // combination of coefficients

      // ignore variables whose coefficient is cancelled by this
      // value of alpha

      if (fabs (ci) < MIN_DENOM)
        continue;

      double
        clbi = clb [indVar],
        cubi = cub [indVar],

        newL = - COUENNE_INFINITY,
        newU =   COUENNE_INFINITY,

        subMin1 = 0., subMin2 = 0.,
        subMax1 = 0., subMax2 = 0.;

#ifdef DEBUG
      printf ("  now tightening x%d\n    denominator: %g * %g + %g * %g = %g\n", 
              (j < n1) ? ind1 [j] : ind2 [j-n1],
              alpha, sa1i, 1-alpha, sa2i, ci);
#endif

      // m[in,ax]Sum[1,2]A must be updated by subtracting the term
      // that is now at the denominator (no need to multiply by
      // signA, sa2 is already adjusted)

      if (clbi <= - COUENNE_INFINITY/10) {

        // we are deleting an infinite bound from the numerator, and
        // whether it goes towards max- or min- Sum? depends on each
        // coefficient

        if (ci > 0.) tickMin++; 
        else         tickMax++;

      } else {

        if (ci > 0.) {

          subMin1 = sa1i * clbi; 
          subMin2 = sa2i * clbi; 

        } else {

          subMax1 = sa1i * clbi;
          subMax2 = sa2i * clbi;
        }
      }

      if (cubi >=   COUENNE_INFINITY/10) {

        // we are deleting an infinite bound from the numerator, and
        // whether it goes towards max- or min- Sum?A it depends on
        // each coefficient

        if (ci > 0.) tickMax++; 
        else         tickMin++;

      } else {

        if (ci > 0.) {

          subMax1 = sa1i * cubi; 
          subMax2 = sa2i * cubi; 

        } else {

          subMin1 = sa1i * cubi;
          subMin2 = sa2i * cubi;
        }
      }

      if (fabs (ci) < MIN_DENOM) { // must check that the numerator behaves

      } else {

        if ((l1 > -COUENNE_INFINITY/10) && 
            (l2 > -COUENNE_INFINITY/10) &&
            (pInfs == tickMax)) {

          newL = ((l1 - l2 - (maxSum1 - maxSum2) + (subMax1 - subMax2)) * alpha + l2 - maxSum2 + subMax2) / ci;

#ifdef DEBUG
          printf ("    \
attempting newL = ((l1 - l2 - (maxSum1 - maxSum2) + (subMax1 - subMax2)) * alpha + l2 - maxSum2 + subMax2) / ci\n    \
((%g - %g - (%g - %g) + (%g - %g)) * %g + %g - %g + %g) / %g = %g\n",
                  l1, l2, maxSum1, maxSum2, subMax1, subMax2, alpha, l2, maxSum2, subMax2, ci, newL);
#endif
        }

        if ((u1 <  COUENNE_INFINITY/10) && 
            (u2 <  COUENNE_INFINITY/10) &&
            (mInfs == tickMin)) {

          newU = ((u1 - u2 - (minSum1 - minSum2) + (subMin1 - subMin2)) * alpha + u2 - minSum2 + subMin2) / ci;

#ifdef DEBUG
          printf ("    \
attempting newU = ((u1 - u2 - (minSum1 - minSum2) + (subMin1 - subMin2)) * alpha + u2 - minSum2 + subMin2) / ci\n    \
((%g - %g - (%g - %g) + (%g - %g)) * %g + %g - %g + %g) / %g = %g\n",
                  u1, u2, minSum1, minSum2, subMin1, subMin2, alpha, u2, minSum2, subMin2, ci, newU);
#endif
        }

        if (ci < 0.) { // should have done the opposite assignment -- just swap them

#ifdef DEBUG
          printf ("    swap'em: %g, %g\n", newL, newU);
#endif

          register double tmp = newL <= - COUENNE_INFINITY / 10 ?   COUENNE_INFINITY : newL;
          newL                = newU >=   COUENNE_INFINITY / 10 ? - COUENNE_INFINITY : newU;
          newU                = tmp;
        }

#ifdef DEBUG
        printf ("    final: %g, %g\n", newL, newU);
#endif
      }

#ifdef DEBUG
      printf ("    bounds for x_%d: [%g,%g] vs [%g,%g]\n", indVar, newL, newU, clb [indVar], cub [indVar]);
#endif

      if ((newL > cubi + COUENNE_EPS) || 
          (newU < clbi - COUENNE_EPS)) {

        delete [] inalphas;
        delete [] alphas;
        delete [] ind1;
        delete [] ind2;

#ifdef DEBUG
        printf ("infeasible\n");
#endif
        return -1;
      }

#ifdef DEBUG
      if (p -> bestSol () &&
          ((p -> bestSol () [indVar] > newU + COUENNE_EPS) ||
           (p -> bestSol () [indVar] < newL - COUENNE_EPS)) &&
          (p -> bestSol () [indVar] >= clbi) &&
          (p -> bestSol () [indVar] <= cubi))

        printf ("optimum violated: %g not in [%g,%g]\n", p -> bestSol () [indVar], newL, newU);
#endif

      double 
        &lbi = Lb [indVar],
        &ubi = Ub [indVar];

      if ((newL > lbi + COUENNE_EPS) && (newL > -COUENNE_INFINITY / 10)) {
                   
        ntightened++;
        lbi = newL;
        newLB. push_back (std::pair <int, double> (indVar, newL));

#ifdef DEBUG
        printf ("    new lower bound: x%d >= %g [%g]\n", indVar, newL, clb [indVar]);
#endif
      }

      if ((newU < ubi - COUENNE_EPS) && (newU < COUENNE_INFINITY / 10)) {

        ntightened++;
        ubi = newU;
        newUB. push_back (std::pair <int, double> (indVar, newU));

#ifdef DEBUG
        printf ("    new upper bound: x%d <= %g [%g]\n", indVar, newU, cub [indVar]);
#endif
      }
    }
  }

#ifdef DEBUG
  printf ("===================================\n");
#endif

  for (std::vector <std::pair <int, double> >::iterator i = newLB. begin (); i != newLB. end (); ++i) {

    if (isInteger [i -> first]) i -> second = ceil (i -> second - COUENNE_EPS);

    if (i -> second > clb [i -> first]) 
      Lb [i -> first] = clb [i -> first] = i -> second;
  }

  for (std::vector <std::pair <int, double> >::iterator i = newUB. begin (); i != newUB. end (); ++i) {

    if (isInteger [i -> first]) i -> second = floor (i -> second + COUENNE_EPS);

    if (i -> second < cub [i -> first]) 
      Ub [i -> first] = cub [i -> first] = i -> second;
  }

  delete [] inalphas;
  delete [] alphas;

  delete [] ind1;
  delete [] ind2;

  return ntightened;
}

}