/* glpios07.c (mixed cover cut generator) */

/***********************************************************************

*  This code is part of GLPK (GNU Linear Programming Kit).

*

*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,

*  2009, 2010 Andrew Makhorin, Department for Applied Informatics,

*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.

*  E-mail: <mao@gnu.org>.

*

*  GLPK is free software: you can redistribute it and/or modify it

*  under the terms of the GNU General Public License as published by

*  the Free Software Foundation, either version 3 of the License, or

*  (at your option) any later version.

*

*  GLPK is distributed in the hope that it will be useful, but WITHOUT

*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

*  License for more details.

*

*  You should have received a copy of the GNU General Public License

*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.

***********************************************************************/

#include "glpios.h"

/*----------------------------------------------------------------------

-- COVER INEQUALITIES

--

-- Consider the set of feasible solutions to 0-1 knapsack problem:

--

--    sum a[j]*x[j] <= b,                                            (1)

--  j in J

--

--    x[j] is binary,                                                (2)

--

-- where, wlog, we assume that a[j] > 0 (since 0-1 variables can be

-- complemented) and a[j] <= b (since a[j] > b implies x[j] = 0).

--

-- A set C within J is called a cover if

--

--    sum a[j] > b.                                                  (3)

--  j in C

--

-- For any cover C the inequality

--

--    sum x[j] <= |C| - 1                                            (4)

--  j in C

--

-- is called a cover inequality and is valid for (1)-(2).

--

-- MIXED COVER INEQUALITIES

--

-- Consider the set of feasible solutions to mixed knapsack problem:

--

--    sum a[j]*x[j] + y <= b,                                        (5)

--  j in J

--

--    x[j] is binary,                                                (6)

--

--    0 <= y <= u is continuous,                                     (7)

--

-- where again we assume that a[j] > 0.

--

-- Let C within J be some set. From (1)-(4) it follows that

--

--    sum a[j] > b - y                                               (8)

--  j in C

--

-- implies

--

--    sum x[j] <= |C| - 1.                                           (9)

--  j in C

--

-- Thus, we need to modify the inequality (9) in such a way that it be

-- a constraint only if the condition (8) is satisfied.

--

-- Consider the following inequality:

--

--    sum x[j] <= |C| - t.                                          (10)

--  j in C

--

-- If 0 < t <= 1, then (10) is equivalent to (9), because all x[j] are

-- binary variables. On the other hand, if t <= 0, (10) being satisfied

-- for any values of x[j] is not a constraint.

--

-- Let

--

--    t' = sum a[j] + y - b.                                        (11)

--       j in C

--

-- It is understood that the condition t' > 0 is equivalent to (8).

-- Besides, from (6)-(7) it follows that t' has an implied upper bound:

--

--    t'max = sum a[j] + u - b.                                     (12)

--          j in C

--

-- This allows to express the parameter t having desired properties:

--

--    t = t' / t'max.                                               (13)

--

-- In fact, t <= 1 by definition, and t > 0 being equivalent to t' > 0

-- is equivalent to (8).

--

-- Thus, the inequality (10), where t is given by formula (13) is valid

-- for (5)-(7).

--

-- Note that if u = 0, then y = 0, so t = 1, and the conditions (8) and

-- (10) is transformed to the conditions (3) and (4).

--

-- GENERATING MIXED COVER CUTS

--

-- To generate a mixed cover cut in the form (10) we need to find such

-- set C which satisfies to the inequality (8) and for which, in turn,

-- the inequality (10) is violated in the current point.

--

-- Substituting t from (13) to (10) gives:

--

--                        1

--    sum x[j] <= |C| - -----  (sum a[j] + y - b),                  (14)

--  j in C              t'max j in C

--

-- and finally we have the cut inequality in the standard form:

--

--    sum x[j] + alfa * y <= beta,                                  (15)

--  j in C

--

-- where:

--

--    alfa = 1 / t'max,                                             (16)

--

--    beta = |C| - alfa *  (sum a[j] - b).                          (17)

--                        j in C                                      */

#if 1

#define MAXTRY 1000

#else

#define MAXTRY 10000

#endif

static int cover2(int n, double a[], double b, double u, double x[],

      double y, int cov[], double *_alfa, double *_beta)

{     /* try to generate mixed cover cut using two-element cover */

      int i, j, try = 0, ret = 0;

      double eps, alfa, beta, temp, rmax = 0.001;

      eps = 0.001 * (1.0 + fabs(b));

      for (i = 0+1; i <= n; i++)

      for (j = i+1; j <= n; j++)

      {  /* C = {i, j} */

         try++;

         if (try > MAXTRY) goto done;

         /* check if condition (8) is satisfied */

         if (a[i] + a[j] + y > b + eps)

         {  /* compute parameters for inequality (15) */

            temp = a[i] + a[j] - b;

            alfa = 1.0 / (temp + u);

            beta = 2.0 - alfa * temp;

            /* compute violation of inequality (15) */

            temp = x[i] + x[j] + alfa * y - beta;

            /* choose C providing maximum violation */

            if (rmax < temp)

            {  rmax = temp;

               cov[1] = i;

               cov[2] = j;

               *_alfa = alfa;

               *_beta = beta;

               ret = 1;

            }

         }

      }

done: return ret;

}

static int cover3(int n, double a[], double b, double u, double x[],

      double y, int cov[], double *_alfa, double *_beta)

{     /* try to generate mixed cover cut using three-element cover */

      int i, j, k, try = 0, ret = 0;

      double eps, alfa, beta, temp, rmax = 0.001;

      eps = 0.001 * (1.0 + fabs(b));

      for (i = 0+1; i <= n; i++)

      for (j = i+1; j <= n; j++)

      for (k = j+1; k <= n; k++)

      {  /* C = {i, j, k} */

         try++;

         if (try > MAXTRY) goto done;

         /* check if condition (8) is satisfied */

         if (a[i] + a[j] + a[k] + y > b + eps)

         {  /* compute parameters for inequality (15) */

            temp = a[i] + a[j] + a[k] - b;

            alfa = 1.0 / (temp + u);

            beta = 3.0 - alfa * temp;

            /* compute violation of inequality (15) */

            temp = x[i] + x[j] + x[k] + alfa * y - beta;

            /* choose C providing maximum violation */

            if (rmax < temp)

            {  rmax = temp;

               cov[1] = i;

               cov[2] = j;

               cov[3] = k;

               *_alfa = alfa;

               *_beta = beta;

               ret = 1;

            }

         }

      }

done: return ret;

}

static int cover4(int n, double a[], double b, double u, double x[],

      double y, int cov[], double *_alfa, double *_beta)

{     /* try to generate mixed cover cut using four-element cover */

      int i, j, k, l, try = 0, ret = 0;

      double eps, alfa, beta, temp, rmax = 0.001;

      eps = 0.001 * (1.0 + fabs(b));

      for (i = 0+1; i <= n; i++)

      for (j = i+1; j <= n; j++)

      for (k = j+1; k <= n; k++)

      for (l = k+1; l <= n; l++)

      {  /* C = {i, j, k, l} */

         try++;

         if (try > MAXTRY) goto done;

         /* check if condition (8) is satisfied */

         if (a[i] + a[j] + a[k] + a[l] + y > b + eps)

         {  /* compute parameters for inequality (15) */

            temp = a[i] + a[j] + a[k] + a[l] - b;

            alfa = 1.0 / (temp + u);

            beta = 4.0 - alfa * temp;

            /* compute violation of inequality (15) */

            temp = x[i] + x[j] + x[k] + x[l] + alfa * y - beta;

            /* choose C providing maximum violation */

            if (rmax < temp)

            {  rmax = temp;

               cov[1] = i;

               cov[2] = j;

               cov[3] = k;

               cov[4] = l;

               *_alfa = alfa;

               *_beta = beta;

               ret = 1;

            }

         }

      }

done: return ret;

}

static int cover(int n, double a[], double b, double u, double x[],

      double y, int cov[], double *alfa, double *beta)

{     /* try to generate mixed cover cut;

         input (see (5)):

         n        is the number of binary variables;

         a[1:n]   are coefficients at binary variables;

         b        is the right-hand side;

         u        is upper bound of continuous variable;

         x[1:n]   are values of binary variables at current point;

         y        is value of continuous variable at current point;

         output (see (15), (16), (17)):

         cov[1:r] are indices of binary variables included in cover C,

                  where r is the set cardinality returned on exit;

         alfa     coefficient at continuous variable;

         beta     is the right-hand side; */

      int j;

      /* perform some sanity checks */

      xassert(n >= 2);

      for (j = 1; j <= n; j++) xassert(a[j] > 0.0);

#if 1 /* ??? */

      xassert(b > -1e-5);

#else

      xassert(b > 0.0);

#endif

      xassert(u >= 0.0);

      for (j = 1; j <= n; j++) xassert(0.0 <= x[j] && x[j] <= 1.0);

      xassert(0.0 <= y && y <= u);

      /* try to generate mixed cover cut */

      if (cover2(n, a, b, u, x, y, cov, alfa, beta)) return 2;

      if (cover3(n, a, b, u, x, y, cov, alfa, beta)) return 3;

      if (cover4(n, a, b, u, x, y, cov, alfa, beta)) return 4;

      return 0;

}

/*----------------------------------------------------------------------

-- lpx_cover_cut - generate mixed cover cut.

--

-- SYNOPSIS

--

-- #include "glplpx.h"

-- int lpx_cover_cut(LPX *lp, int len, int ind[], double val[],

--    double work[]);

--

-- DESCRIPTION

--

-- The routine lpx_cover_cut generates a mixed cover cut for a given

-- row of the MIP problem.

--

-- The given row of the MIP problem should be explicitly specified in

-- the form:

--

--    sum{j in J} a[j]*x[j] <= b.                                    (1)

--

-- On entry indices (ordinal numbers) of structural variables, which

-- have non-zero constraint coefficients, should be placed in locations

-- ind[1], ..., ind[len], and corresponding constraint coefficients

-- should be placed in locations val[1], ..., val[len]. The right-hand

-- side b should be stored in location val[0].

--

-- The working array work should have at least nb locations, where nb

-- is the number of binary variables in (1).

--

-- The routine generates a mixed cover cut in the same form as (1) and

-- stores the cut coefficients and right-hand side in the same way as

-- just described above.

--

-- RETURNS

--

-- If the cutting plane has been successfully generated, the routine

-- returns 1 <= len' <= n, which is the number of non-zero coefficients

-- in the inequality constraint. Otherwise, the routine returns zero. */

static int lpx_cover_cut(LPX *lp, int len, int ind[], double val[],

      double work[])

{     int cov[1+4], j, k, nb, newlen, r;

      double f_min, f_max, alfa, beta, u, *x = work, y;

      /* substitute and remove fixed variables */

      newlen = 0;

      for (k = 1; k <= len; k++)

      {  j = ind[k];

         if (lpx_get_col_type(lp, j) == LPX_FX)

            val[0] -= val[k] * lpx_get_col_lb(lp, j);

         else

         {  newlen++;

            ind[newlen] = ind[k];

            val[newlen] = val[k];

         }

      }

      len = newlen;

      /* move binary variables to the beginning of the list so that

         elements 1, 2, ..., nb correspond to binary variables, and

         elements nb+1, nb+2, ..., len correspond to rest variables */

      nb = 0;

      for (k = 1; k <= len; k++)

      {  j = ind[k];

         if (lpx_get_col_kind(lp, j) == LPX_IV &&

             lpx_get_col_type(lp, j) == LPX_DB &&

             lpx_get_col_lb(lp, j) == 0.0 &&

             lpx_get_col_ub(lp, j) == 1.0)

         {  /* binary variable */

            int ind_k;

            double val_k;

            nb++;

            ind_k = ind[nb], val_k = val[nb];

            ind[nb] = ind[k], val[nb] = val[k];

            ind[k] = ind_k, val[k] = val_k;

         }

      }

      /* now the specified row has the form:

         sum a[j]*x[j] + sum a[j]*y[j] <= b,

         where x[j] are binary variables, y[j] are rest variables */

      /* at least two binary variables are needed */

      if (nb < 2) return 0;

      /* compute implied lower and upper bounds for sum a[j]*y[j] */

      f_min = f_max = 0.0;

      for (k = nb+1; k <= len; k++)

      {  j = ind[k];

         /* both bounds must be finite */

         if (lpx_get_col_type(lp, j) != LPX_DB) return 0;

         if (val[k] > 0.0)

         {  f_min += val[k] * lpx_get_col_lb(lp, j);

            f_max += val[k] * lpx_get_col_ub(lp, j);

         }

         else

         {  f_min += val[k] * lpx_get_col_ub(lp, j);

            f_max += val[k] * lpx_get_col_lb(lp, j);

         }

      }

      /* sum a[j]*x[j] + sum a[j]*y[j] <= b ===>

         sum a[j]*x[j] + (sum a[j]*y[j] - f_min) <= b - f_min ===>

         sum a[j]*x[j] + y <= b - f_min,

         where y = sum a[j]*y[j] - f_min;

         note that 0 <= y <= u, u = f_max - f_min */

      /* determine upper bound of y */

      u = f_max - f_min;

      /* determine value of y at the current point */

      y = 0.0;

      for (k = nb+1; k <= len; k++)

      {  j = ind[k];

         y += val[k] * lpx_get_col_prim(lp, j);

      }

      y -= f_min;

      if (y < 0.0) y = 0.0;

      if (y > u) y = u;

      /* modify the right-hand side b */

      val[0] -= f_min;

      /* now the transformed row has the form:

         sum a[j]*x[j] + y <= b, where 0 <= y <= u */

      /* determine values of x[j] at the current point */

      for (k = 1; k <= nb; k++)

      {  j = ind[k];

         x[k] = lpx_get_col_prim(lp, j);

         if (x[k] < 0.0) x[k] = 0.0;

         if (x[k] > 1.0) x[k] = 1.0;

      }

      /* if a[j] < 0, replace x[j] by its complement 1 - x'[j] */

      for (k = 1; k <= nb; k++)

      {  if (val[k] < 0.0)

         {  ind[k] = - ind[k];

            val[k] = - val[k];

            val[0] += val[k];

            x[k] = 1.0 - x[k];

         }

      }

      /* try to generate a mixed cover cut for the transformed row */

      r = cover(nb, val, val[0], u, x, y, cov, &alfa, &beta);

      if (r == 0) return 0;

      xassert(2 <= r && r <= 4);

      /* now the cut is in the form:

         sum{j in C} x[j] + alfa * y <= beta */

      /* store the right-hand side beta */

      ind[0] = 0, val[0] = beta;

      /* restore the original ordinal numbers of x[j] */

      for (j = 1; j <= r; j++) cov[j] = ind[cov[j]];

      /* store cut coefficients at binary variables complementing back

         the variables having negative row coefficients */

      xassert(r <= nb);

      for (k = 1; k <= r; k++)

      {  if (cov[k] > 0)

         {  ind[k] = +cov[k];

            val[k] = +1.0;

         }

         else

         {  ind[k] = -cov[k];

            val[k] = -1.0;

            val[0] -= 1.0;

         }

      }

      /* substitute y = sum a[j]*y[j] - f_min */

      for (k = nb+1; k <= len; k++)

      {  r++;

         ind[r] = ind[k];

         val[r] = alfa * val[k];

      }

      val[0] += alfa * f_min;

      xassert(r <= len);

      len = r;

      return len;

}

/*----------------------------------------------------------------------

-- lpx_eval_row - compute explictily specified row.

--

-- SYNOPSIS

--

-- #include "glplpx.h"

-- double lpx_eval_row(LPX *lp, int len, int ind[], double val[]);

--

-- DESCRIPTION

--

-- The routine lpx_eval_row computes the primal value of an explicitly

-- specified row using current values of structural variables.

--

-- The explicitly specified row may be thought as a linear form:

--

--    y = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n],

--

-- where y is an auxiliary variable for this row, a[j] are coefficients

-- of the linear form, x[m+j] are structural variables.

--

-- On entry column indices and numerical values of non-zero elements of

-- the row should be stored in locations ind[1], ..., ind[len] and

-- val[1], ..., val[len], where len is the number of non-zero elements.

-- The array ind and val are not changed on exit.

--

-- RETURNS

--

-- The routine returns a computed value of y, the auxiliary variable of

-- the specified row. */

static double lpx_eval_row(LPX *lp, int len, int ind[], double val[])

{     int n = lpx_get_num_cols(lp);

      int j, k;

      double sum = 0.0;

      if (len < 0)

         xerror("lpx_eval_row: len = %d; invalid row length\n", len);

      for (k = 1; k <= len; k++)

      {  j = ind[k];

         if (!(1 <= j && j <= n))

            xerror("lpx_eval_row: j = %d; column number out of range\n",

               j);

         sum += val[k] * lpx_get_col_prim(lp, j);

      }

      return sum;

}

/***********************************************************************

*  NAME

*

*  ios_cov_gen - generate mixed cover cuts

*

*  SYNOPSIS

*

*  #include "glpios.h"

*  void ios_cov_gen(glp_tree *tree);

*

*  DESCRIPTION

*

*  The routine ios_cov_gen generates mixed cover cuts for the current

*  point and adds them to the cut pool. */

void ios_cov_gen(glp_tree *tree)

{     glp_prob *prob = tree->mip;

      int m = lpx_get_num_rows(prob);

      int n = lpx_get_num_cols(prob);

      int i, k, type, kase, len, *ind;

      double r, *val, *work;

      xassert(lpx_get_status(prob) == LPX_OPT);

      /* allocate working arrays */

      ind = xcalloc(1+n, sizeof(int));

      val = xcalloc(1+n, sizeof(double));

      work = xcalloc(1+n, sizeof(double));

      /* look through all rows */

      for (i = 1; i <= m; i++)

      for (kase = 1; kase <= 2; kase++)

      {  type = lpx_get_row_type(prob, i);

         if (kase == 1)

         {  /* consider rows of '<=' type */

            if (!(type == LPX_UP || type == LPX_DB)) continue;

            len = lpx_get_mat_row(prob, i, ind, val);

            val[0] = lpx_get_row_ub(prob, i);

         }

         else

         {  /* consider rows of '>=' type */

            if (!(type == LPX_LO || type == LPX_DB)) continue;

            len = lpx_get_mat_row(prob, i, ind, val);

            for (k = 1; k <= len; k++) val[k] = - val[k];

            val[0] = - lpx_get_row_lb(prob, i);

         }

         /* generate mixed cover cut:

            sum{j in J} a[j] * x[j] <= b */

         len = lpx_cover_cut(prob, len, ind, val, work);

         if (len == 0) continue;

         /* at the current point the cut inequality is violated, i.e.

            sum{j in J} a[j] * x[j] - b > 0 */

         r = lpx_eval_row(prob, len, ind, val) - val[0];

         if (r < 1e-3) continue;

         /* add the cut to the cut pool */

         glp_ios_add_row(tree, NULL, GLP_RF_COV, 0, len, ind, val,

            GLP_UP, val[0]);

      }

      /* free working arrays */

      xfree(ind);

      xfree(val);

      xfree(work);

      return;

}

/* eof */