/* glpios06.c (MIR 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"

#define _MIR_DEBUG 0

#define MAXAGGR 5

/* maximal number of rows which can be aggregated */

struct MIR

{     /* MIR cut generator working area */

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

      /* global information valid for the root subproblem */

      int m;

      /* number of rows (in the root subproblem) */

      int n;

      /* number of columns */

      char *skip; /* char skip[1+m]; */

      /* skip[i], 1 <= i <= m, is a flag that means that row i should

         not be used because (1) it is not suitable, or (2) because it

         has been used in the aggregated constraint */

      char *isint; /* char isint[1+m+n]; */

      /* isint[k], 1 <= k <= m+n, is a flag that means that variable

         x[k] is integer (otherwise, continuous) */

      double *lb; /* double lb[1+m+n]; */

      /* lb[k], 1 <= k <= m+n, is lower bound of x[k]; -DBL_MAX means

         that x[k] has no lower bound */

      int *vlb; /* int vlb[1+m+n]; */

      /* vlb[k] = k', 1 <= k <= m+n, is the number of integer variable,

         which defines variable lower bound x[k] >= lb[k] * x[k']; zero

         means that x[k] has simple lower bound */

      double *ub; /* double ub[1+m+n]; */

      /* ub[k], 1 <= k <= m+n, is upper bound of x[k]; +DBL_MAX means

         that x[k] has no upper bound */

      int *vub; /* int vub[1+m+n]; */

      /* vub[k] = k', 1 <= k <= m+n, is the number of integer variable,

         which defines variable upper bound x[k] <= ub[k] * x[k']; zero

         means that x[k] has simple upper bound */

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

      /* current (fractional) point to be separated */

      double *x; /* double x[1+m+n]; */

      /* x[k] is current value of auxiliary (1 <= k <= m) or structural

         (m+1 <= k <= m+n) variable */

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

      /* aggregated constraint sum a[k] * x[k] = b, which is a linear

         combination of original constraints transformed to equalities

         by introducing auxiliary variables */

      int agg_cnt;

      /* number of rows (original constraints) used to build aggregated

         constraint, 1 <= agg_cnt <= MAXAGGR */

      int *agg_row; /* int agg_row[1+MAXAGGR]; */

      /* agg_row[k], 1 <= k <= agg_cnt, is the row number used to build

         aggregated constraint */

      IOSVEC *agg_vec; /* IOSVEC agg_vec[1:m+n]; */

      /* sparse vector of aggregated constraint coefficients, a[k] */

      double agg_rhs;

      /* right-hand side of the aggregated constraint, b */

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

      /* bound substitution flags for modified constraint */

      char *subst; /* char subst[1+m+n]; */

      /* subst[k], 1 <= k <= m+n, is a bound substitution flag used for

         variable x[k]:

         '?' - x[k] is missing in modified constraint

         'L' - x[k] = (lower bound) + x'[k]

         'U' - x[k] = (upper bound) - x'[k] */

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

      /* modified constraint sum a'[k] * x'[k] = b', where x'[k] >= 0,

         derived from aggregated constraint by substituting bounds;

         note that due to substitution of variable bounds there may be

         additional terms in the modified constraint */

      IOSVEC *mod_vec; /* IOSVEC mod_vec[1:m+n]; */

      /* sparse vector of modified constraint coefficients, a'[k] */

      double mod_rhs;

      /* right-hand side of the modified constraint, b' */

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

      /* cutting plane sum alpha[k] * x[k] <= beta */

      IOSVEC *cut_vec; /* IOSVEC cut_vec[1:m+n]; */

      /* sparse vector of cutting plane coefficients, alpha[k] */

      double cut_rhs;

      /* right-hand size of the cutting plane, beta */

};

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

*  NAME

*

*  ios_mir_init - initialize MIR cut generator

*

*  SYNOPSIS

*

*  #include "glpios.h"

*  void *ios_mir_init(glp_tree *tree);

*

*  DESCRIPTION

*

*  The routine ios_mir_init initializes the MIR cut generator assuming

*  that the current subproblem is the root subproblem.

*

*  RETURNS

*

*  The routine ios_mir_init returns a pointer to the MIR cut generator

*  working area. */

static void set_row_attrib(glp_tree *tree, struct MIR *mir)

{     /* set global row attributes */

      glp_prob *mip = tree->mip;

      int m = mir->m;

      int k;

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

      {  GLPROW *row = mip->row[k];

         mir->skip[k] = 0;

         mir->isint[k] = 0;

         switch (row->type)

         {  case GLP_FR:

               mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;

            case GLP_LO:

               mir->lb[k] = row->lb, mir->ub[k] = +DBL_MAX; break;

            case GLP_UP:

               mir->lb[k] = -DBL_MAX, mir->ub[k] = row->ub; break;

            case GLP_DB:

               mir->lb[k] = row->lb, mir->ub[k] = row->ub; break;

            case GLP_FX:

               mir->lb[k] = mir->ub[k] = row->lb; break;

            default:

               xassert(row != row);

         }

         mir->vlb[k] = mir->vub[k] = 0;

      }

      return;

}

static void set_col_attrib(glp_tree *tree, struct MIR *mir)

{     /* set global column attributes */

      glp_prob *mip = tree->mip;

      int m = mir->m;

      int n = mir->n;

      int k;

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

      {  GLPCOL *col = mip->col[k-m];

         switch (col->kind)

         {  case GLP_CV:

               mir->isint[k] = 0; break;

            case GLP_IV:

               mir->isint[k] = 1; break;

            default:

               xassert(col != col);

         }

         switch (col->type)

         {  case GLP_FR:

               mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;

            case GLP_LO:

               mir->lb[k] = col->lb, mir->ub[k] = +DBL_MAX; break;

            case GLP_UP:

               mir->lb[k] = -DBL_MAX, mir->ub[k] = col->ub; break;

            case GLP_DB:

               mir->lb[k] = col->lb, mir->ub[k] = col->ub; break;

            case GLP_FX:

               mir->lb[k] = mir->ub[k] = col->lb; break;

            default:

               xassert(col != col);

         }

         mir->vlb[k] = mir->vub[k] = 0;

      }

      return;

}

static void set_var_bounds(glp_tree *tree, struct MIR *mir)

{     /* set variable bounds */

      glp_prob *mip = tree->mip;

      int m = mir->m;

      GLPAIJ *aij;

      int i, k1, k2;

      double a1, a2;

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

      {  /* we need the row to be '>= 0' or '<= 0' */

         if (!(mir->lb[i] == 0.0 && mir->ub[i] == +DBL_MAX ||

               mir->lb[i] == -DBL_MAX && mir->ub[i] == 0.0)) continue;

         /* take first term */

         aij = mip->row[i]->ptr;

         if (aij == NULL) continue;

         k1 = m + aij->col->j, a1 = aij->val;

         /* take second term */

         aij = aij->r_next;

         if (aij == NULL) continue;

         k2 = m + aij->col->j, a2 = aij->val;

         /* there must be only two terms */

         if (aij->r_next != NULL) continue;

         /* interchange terms, if needed */

         if (!mir->isint[k1] && mir->isint[k2])

            ;

         else if (mir->isint[k1] && !mir->isint[k2])

         {  k2 = k1, a2 = a1;

            k1 = m + aij->col->j, a1 = aij->val;

         }

         else

         {  /* both terms are either continuous or integer */

            continue;

         }

         /* x[k2] should be double-bounded */

         if (mir->lb[k2] == -DBL_MAX || mir->ub[k2] == +DBL_MAX ||

             mir->lb[k2] == mir->ub[k2]) continue;

         /* change signs, if necessary */

         if (mir->ub[i] == 0.0) a1 = - a1, a2 = - a2;

         /* now the row has the form a1 * x1 + a2 * x2 >= 0, where x1

            is continuous, x2 is integer */

         if (a1 > 0.0)

         {  /* x1 >= - (a2 / a1) * x2 */

            if (mir->vlb[k1] == 0)

            {  /* set variable lower bound for x1 */

               mir->lb[k1] = - a2 / a1;

               mir->vlb[k1] = k2;

               /* the row should not be used */

               mir->skip[i] = 1;

            }

         }

         else /* a1 < 0.0 */

         {  /* x1 <= - (a2 / a1) * x2 */

            if (mir->vub[k1] == 0)

            {  /* set variable upper bound for x1 */

               mir->ub[k1] = - a2 / a1;

               mir->vub[k1] = k2;

               /* the row should not be used */

               mir->skip[i] = 1;

            }

         }

      }

      return;

}

static void mark_useless_rows(glp_tree *tree, struct MIR *mir)

{     /* mark rows which should not be used */

      glp_prob *mip = tree->mip;

      int m = mir->m;

      GLPAIJ *aij;

      int i, k, nv;

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

      {  /* free rows should not be used */

         if (mir->lb[i] == -DBL_MAX && mir->ub[i] == +DBL_MAX)

         {  mir->skip[i] = 1;

            continue;

         }

         nv = 0;

         for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)

         {  k = m + aij->col->j;

            /* rows with free variables should not be used */

            if (mir->lb[k] == -DBL_MAX && mir->ub[k] == +DBL_MAX)

            {  mir->skip[i] = 1;

               break;

            }

            /* rows with integer variables having infinite (lower or

               upper) bound should not be used */

            if (mir->isint[k] && mir->lb[k] == -DBL_MAX ||

                mir->isint[k] && mir->ub[k] == +DBL_MAX)

            {  mir->skip[i] = 1;

               break;

            }

            /* count non-fixed variables */

            if (!(mir->vlb[k] == 0 && mir->vub[k] == 0 &&

                  mir->lb[k] == mir->ub[k])) nv++;

         }

         /* rows with all variables fixed should not be used */

         if (nv == 0)

         {  mir->skip[i] = 1;

            continue;

         }

      }

      return;

}

void *ios_mir_init(glp_tree *tree)

{     /* initialize MIR cut generator */

      glp_prob *mip = tree->mip;

      int m = mip->m;

      int n = mip->n;

      struct MIR *mir;

#if _MIR_DEBUG

      xprintf("ios_mir_init: warning: debug mode enabled\n");

#endif

      /* allocate working area */

      mir = xmalloc(sizeof(struct MIR));

      mir->m = m;

      mir->n = n;

      mir->skip = xcalloc(1+m, sizeof(char));

      mir->isint = xcalloc(1+m+n, sizeof(char));

      mir->lb = xcalloc(1+m+n, sizeof(double));

      mir->vlb = xcalloc(1+m+n, sizeof(int));

      mir->ub = xcalloc(1+m+n, sizeof(double));

      mir->vub = xcalloc(1+m+n, sizeof(int));

      mir->x = xcalloc(1+m+n, sizeof(double));

      mir->agg_row = xcalloc(1+MAXAGGR, sizeof(int));

      mir->agg_vec = ios_create_vec(m+n);

      mir->subst = xcalloc(1+m+n, sizeof(char));

      mir->mod_vec = ios_create_vec(m+n);

      mir->cut_vec = ios_create_vec(m+n);

      /* set global row attributes */

      set_row_attrib(tree, mir);

      /* set global column attributes */

      set_col_attrib(tree, mir);

      /* set variable bounds */

      set_var_bounds(tree, mir);

      /* mark rows which should not be used */

      mark_useless_rows(tree, mir);

      return mir;

}

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

*  NAME

*

*  ios_mir_gen - generate MIR cuts

*

*  SYNOPSIS

*

*  #include "glpios.h"

*  void ios_mir_gen(glp_tree *tree, void *gen, IOSPOOL *pool);

*

*  DESCRIPTION

*

*  The routine ios_mir_gen generates MIR cuts for the current point and

*  adds them to the cut pool. */

static void get_current_point(glp_tree *tree, struct MIR *mir)

{     /* obtain current point */

      glp_prob *mip = tree->mip;

      int m = mir->m;

      int n = mir->n;

      int k;

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

         mir->x[k] = mip->row[k]->prim;

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

         mir->x[k] = mip->col[k-m]->prim;

      return;

}

#if _MIR_DEBUG

static void check_current_point(struct MIR *mir)

{     /* check current point */

      int m = mir->m;

      int n = mir->n;

      int k, kk;

      double lb, ub, eps;

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

      {  /* determine lower bound */

         lb = mir->lb[k];

         kk = mir->vlb[k];

         if (kk != 0)

         {  xassert(lb != -DBL_MAX);

            xassert(!mir->isint[k]);

            xassert(mir->isint[kk]);

            lb *= mir->x[kk];

         }

         /* check lower bound */

         if (lb != -DBL_MAX)

         {  eps = 1e-6 * (1.0 + fabs(lb));

            xassert(mir->x[k] >= lb - eps);

         }

         /* determine upper bound */

         ub = mir->ub[k];

         kk = mir->vub[k];

         if (kk != 0)

         {  xassert(ub != +DBL_MAX);

            xassert(!mir->isint[k]);

            xassert(mir->isint[kk]);

            ub *= mir->x[kk];

         }

         /* check upper bound */

         if (ub != +DBL_MAX)

         {  eps = 1e-6 * (1.0 + fabs(ub));

            xassert(mir->x[k] <= ub + eps);

         }

      }

      return;

}

#endif

static void initial_agg_row(glp_tree *tree, struct MIR *mir, int i)

{     /* use original i-th row as initial aggregated constraint */

      glp_prob *mip = tree->mip;

      int m = mir->m;

      GLPAIJ *aij;

      xassert(1 <= i && i <= m);

      xassert(!mir->skip[i]);

      /* mark i-th row in order not to use it in the same aggregated

         constraint */

      mir->skip[i] = 2;

      mir->agg_cnt = 1;

      mir->agg_row[1] = i;

      /* use x[i] - sum a[i,j] * x[m+j] = 0, where x[i] is auxiliary

         variable of row i, x[m+j] are structural variables */

      ios_clear_vec(mir->agg_vec);

      ios_set_vj(mir->agg_vec, i, 1.0);

      for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)

         ios_set_vj(mir->agg_vec, m + aij->col->j, - aij->val);

      mir->agg_rhs = 0.0;

#if _MIR_DEBUG

      ios_check_vec(mir->agg_vec);

#endif

      return;

}

#if _MIR_DEBUG

static void check_agg_row(struct MIR *mir)

{     /* check aggregated constraint */

      int m = mir->m;

      int n = mir->n;

      int j, k;

      double r, big;

      /* compute the residual r = sum a[k] * x[k] - b and determine

         big = max(1, |a[k]|, |b|) */

      r = 0.0, big = 1.0;

      for (j = 1; j <= mir->agg_vec->nnz; j++)

      {  k = mir->agg_vec->ind[j];

         xassert(1 <= k && k <= m+n);

         r += mir->agg_vec->val[j] * mir->x[k];

         if (big < fabs(mir->agg_vec->val[j]))

            big = fabs(mir->agg_vec->val[j]);

      }

      r -= mir->agg_rhs;

      if (big < fabs(mir->agg_rhs))

         big = fabs(mir->agg_rhs);

      /* the residual must be close to zero */

      xassert(fabs(r) <= 1e-6 * big);

      return;

}

#endif

static void subst_fixed_vars(struct MIR *mir)

{     /* substitute fixed variables into aggregated constraint */

      int m = mir->m;

      int n = mir->n;

      int j, k;

      for (j = 1; j <= mir->agg_vec->nnz; j++)

      {  k = mir->agg_vec->ind[j];

         xassert(1 <= k && k <= m+n);

         if (mir->vlb[k] == 0 && mir->vub[k] == 0 &&

             mir->lb[k] == mir->ub[k])

         {  /* x[k] is fixed */

            mir->agg_rhs -= mir->agg_vec->val[j] * mir->lb[k];

            mir->agg_vec->val[j] = 0.0;

         }

      }

      /* remove terms corresponding to fixed variables */

      ios_clean_vec(mir->agg_vec, DBL_EPSILON);

#if _MIR_DEBUG

      ios_check_vec(mir->agg_vec);

#endif

      return;

}

static void bound_subst_heur(struct MIR *mir)

{     /* bound substitution heuristic */

      int m = mir->m;

      int n = mir->n;

      int j, k, kk;

      double d1, d2;

      for (j = 1; j <= mir->agg_vec->nnz; j++)

      {  k = mir->agg_vec->ind[j];

         xassert(1 <= k && k <= m+n);

         if (mir->isint[k]) continue; /* skip integer variable */

         /* compute distance from x[k] to its lower bound */

         kk = mir->vlb[k];

         if (kk == 0)

         {  if (mir->lb[k] == -DBL_MAX)

               d1 = DBL_MAX;

            else

               d1 = mir->x[k] - mir->lb[k];

         }

         else

         {  xassert(1 <= kk && kk <= m+n);

            xassert(mir->isint[kk]);

            xassert(mir->lb[k] != -DBL_MAX);

            d1 = mir->x[k] - mir->lb[k] * mir->x[kk];

         }

         /* compute distance from x[k] to its upper bound */

         kk = mir->vub[k];

         if (kk == 0)

         {  if (mir->vub[k] == +DBL_MAX)

               d2 = DBL_MAX;

            else

               d2 = mir->ub[k] - mir->x[k];

         }

         else

         {  xassert(1 <= kk && kk <= m+n);

            xassert(mir->isint[kk]);

            xassert(mir->ub[k] != +DBL_MAX);

            d2 = mir->ub[k] * mir->x[kk] - mir->x[k];

         }

         /* x[k] cannot be free */

         xassert(d1 != DBL_MAX || d2 != DBL_MAX);

         /* choose the bound which is closer to x[k] */

         xassert(mir->subst[k] == '?');

         if (d1 <= d2)

            mir->subst[k] = 'L';

         else

            mir->subst[k] = 'U';

      }

      return;

}

static void build_mod_row(struct MIR *mir)

{     /* substitute bounds and build modified constraint */

      int m = mir->m;

      int n = mir->n;

      int j, jj, k, kk;

      /* initially modified constraint is aggregated constraint */

      ios_copy_vec(mir->mod_vec, mir->agg_vec);

      mir->mod_rhs = mir->agg_rhs;

#if _MIR_DEBUG

      ios_check_vec(mir->mod_vec);

#endif

      /* substitute bounds for continuous variables; note that due to

         substitution of variable bounds additional terms may appear in

         modified constraint */

      for (j = mir->mod_vec->nnz; j >= 1; j--)

      {  k = mir->mod_vec->ind[j];

         xassert(1 <= k && k <= m+n);

         if (mir->isint[k]) continue; /* skip integer variable */

         if (mir->subst[k] == 'L')

         {  /* x[k] = (lower bound) + x'[k] */

            xassert(mir->lb[k] != -DBL_MAX);

            kk = mir->vlb[k];

            if (kk == 0)

            {  /* x[k] = lb[k] + x'[k] */

               mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];

            }

            else

            {  /* x[k] = lb[k] * x[kk] + x'[k] */

               xassert(mir->isint[kk]);

               jj = mir->mod_vec->pos[kk];

               if (jj == 0)

               {  ios_set_vj(mir->mod_vec, kk, 1.0);

                  jj = mir->mod_vec->pos[kk];

                  mir->mod_vec->val[jj] = 0.0;

               }

               mir->mod_vec->val[jj] +=

                  mir->mod_vec->val[j] * mir->lb[k];

            }

         }

         else if (mir->subst[k] == 'U')

         {  /* x[k] = (upper bound) - x'[k] */

            xassert(mir->ub[k] != +DBL_MAX);

            kk = mir->vub[k];

            if (kk == 0)

            {  /* x[k] = ub[k] - x'[k] */

               mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];

            }

            else

            {  /* x[k] = ub[k] * x[kk] - x'[k] */

               xassert(mir->isint[kk]);

               jj = mir->mod_vec->pos[kk];

               if (jj == 0)

               {  ios_set_vj(mir->mod_vec, kk, 1.0);

                  jj = mir->mod_vec->pos[kk];

                  mir->mod_vec->val[jj] = 0.0;

               }

               mir->mod_vec->val[jj] +=

                  mir->mod_vec->val[j] * mir->ub[k];

            }

            mir->mod_vec->val[j] = - mir->mod_vec->val[j];

         }

         else

            xassert(k != k);

      }

#if _MIR_DEBUG

      ios_check_vec(mir->mod_vec);

#endif

      /* substitute bounds for integer variables */

      for (j = 1; j <= mir->mod_vec->nnz; j++)

      {  k = mir->mod_vec->ind[j];

         xassert(1 <= k && k <= m+n);

         if (!mir->isint[k]) continue; /* skip continuous variable */

         xassert(mir->subst[k] == '?');

         xassert(mir->vlb[k] == 0 && mir->vub[k] == 0);

         xassert(mir->lb[k] != -DBL_MAX && mir->ub[k] != +DBL_MAX);

         if (fabs(mir->lb[k]) <= fabs(mir->ub[k]))

         {  /* x[k] = lb[k] + x'[k] */

            mir->subst[k] = 'L';

            mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];

         }

         else

         {  /* x[k] = ub[k] - x'[k] */

            mir->subst[k] = 'U';

            mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];

            mir->mod_vec->val[j] = - mir->mod_vec->val[j];

         }

      }

#if _MIR_DEBUG

      ios_check_vec(mir->mod_vec);

#endif

      return;

}

#if _MIR_DEBUG

static void check_mod_row(struct MIR *mir)

{     /* check modified constraint */

      int m = mir->m;

      int n = mir->n;

      int j, k, kk;

      double r, big, x;

      /* compute the residual r = sum a'[k] * x'[k] - b' and determine

         big = max(1, |a[k]|, |b|) */

      r = 0.0, big = 1.0;

      for (j = 1; j <= mir->mod_vec->nnz; j++)

      {  k = mir->mod_vec->ind[j];

         xassert(1 <= k && k <= m+n);

         if (mir->subst[k] == 'L')

         {  /* x'[k] = x[k] - (lower bound) */

            xassert(mir->lb[k] != -DBL_MAX);

            kk = mir->vlb[k];

            if (kk == 0)

               x = mir->x[k] - mir->lb[k];

            else

               x = mir->x[k] - mir->lb[k] * mir->x[kk];

         }

         else if (mir->subst[k] == 'U')

         {  /* x'[k] = (upper bound) - x[k] */

            xassert(mir->ub[k] != +DBL_MAX);

            kk = mir->vub[k];

            if (kk == 0)

               x = mir->ub[k] - mir->x[k];

            else

               x = mir->ub[k] * mir->x[kk] - mir->x[k];

         }

         else

            xassert(k != k);

         r += mir->mod_vec->val[j] * x;

         if (big < fabs(mir->mod_vec->val[j]))

            big = fabs(mir->mod_vec->val[j]);

      }

      r -= mir->mod_rhs;

      if (big < fabs(mir->mod_rhs))

         big = fabs(mir->mod_rhs);

      /* the residual must be close to zero */

      xassert(fabs(r) <= 1e-6 * big);

      return;

}

#endif

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

*  mir_ineq - construct MIR inequality

*

*  Given the single constraint mixed integer set

*

*                    |N|

*     X = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s},

*                    +      +  j in N

*

*  this routine constructs the mixed integer rounding (MIR) inequality

*

*     sum   alpha[j] * x[j] <= beta + gamma * s,

*    j in N

*

*  which is valid for X.

*

*  If the MIR inequality has been successfully constructed, the routine

*  returns zero. Otherwise, if b is close to nearest integer, there may

*  be numeric difficulties due to big coefficients; so in this case the

*  routine returns non-zero. */

static int mir_ineq(const int n, const double a[], const double b,

      double alpha[], double *beta, double *gamma)

{     int j;

      double f, t;

      if (fabs(b - floor(b + .5)) < 0.01)

         return 1;

      f = b - floor(b);

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

      {  t = (a[j] - floor(a[j])) - f;

         if (t <= 0.0)

            alpha[j] = floor(a[j]);

         else

            alpha[j] = floor(a[j]) + t / (1.0 - f);

      }

      *beta = floor(b);

      *gamma = 1.0 / (1.0 - f);

      return 0;

}

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

*  cmir_ineq - construct c-MIR inequality

*

*  Given the mixed knapsack set

*

*      MK              |N|

*     X   = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s,

*                      +      +  j in N

*

*             x[j] <= u[j]},

*

*  a subset C of variables to be complemented, and a divisor delta > 0,

*  this routine constructs the complemented MIR (c-MIR) inequality

*

*     sum   alpha[j] * x[j] <= beta + gamma * s,

*    j in N

*                      MK

*  which is valid for X  .

*

*  If the c-MIR inequality has been successfully constructed, the

*  routine returns zero. Otherwise, if there is a risk of numerical

*  difficulties due to big coefficients (see comments to the routine

*  mir_ineq), the routine cmir_ineq returns non-zero. */

static int cmir_ineq(const int n, const double a[], const double b,

      const double u[], const char cset[], const double delta,

      double alpha[], double *beta, double *gamma)

{     int j;

      double *aa, bb;

      aa = alpha, bb = b;

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

      {  aa[j] = a[j] / delta;

         if (cset[j])

            aa[j] = - aa[j], bb -= a[j] * u[j];

      }

      bb /= delta;

      if (mir_ineq(n, aa, bb, alpha, beta, gamma)) return 1;

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

      {  if (cset[j])

            alpha[j] = - alpha[j], *beta += alpha[j] * u[j];

      }

      *gamma /= delta;

      return 0;

}

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

*  cmir_sep - c-MIR separation heuristic

*

*  Given the mixed knapsack set

*

*      MK              |N|

*     X   = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s,

*                      +      +  j in N

*

*             x[j] <= u[j]}

*

*                           *   *

*  and a fractional point (x , s ), this routine tries to construct

*  c-MIR inequality

*

*     sum   alpha[j] * x[j] <= beta + gamma * s,

*    j in N

*                      MK

*  which is valid for X   and has (desirably maximal) violation at the

*  fractional point given. This is attained by choosing an appropriate

*  set C of variables to be complemented and a divisor delta > 0, which

*  together define corresponding c-MIR inequality.

*

*  If a violated c-MIR inequality has been successfully constructed,

*  the routine returns its violation:

*

*                       *                      *

*     sum   alpha[j] * x [j] - beta - gamma * s ,

*    j in N

*

*  which is positive. In case of failure the routine returns zero. */

struct vset { int j; double v; };

static int cmir_cmp(const void *p1, const void *p2)

{     const struct vset *v1 = p1, *v2 = p2;

      if (v1->v < v2->v) return -1;

      if (v1->v > v2->v) return +1;

      return 0;

}

static double cmir_sep(const int n, const double a[], const double b,

      const double u[], const double x[], const double s,

      double alpha[], double *beta, double *gamma)

{     int fail, j, k, nv, v;

      double delta, eps, d_try[1+3], r, r_best;

      char *cset;

      struct vset *vset;

      /* allocate working arrays */

      cset = xcalloc(1+n, sizeof(char));

      vset = xcalloc(1+n, sizeof(struct vset));

      /* choose initial C */

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

         cset[j] = (char)(x[j] >= 0.5 * u[j]);

      /* choose initial delta */

      r_best = delta = 0.0;

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

      {  xassert(a[j] != 0.0);

         /* if x[j] is close to its bounds, skip it */

         eps = 1e-9 * (1.0 + fabs(u[j]));

         if (x[j] < eps || x[j] > u[j] - eps) continue;

         /* try delta = |a[j]| to construct c-MIR inequality */

         fail = cmir_ineq(n, a, b, u, cset, fabs(a[j]), alpha, beta,

            gamma);

         if (fail) continue;

         /* compute violation */

         r = - (*beta) - (*gamma) * s;

         for (k = 1; k <= n; k++) r += alpha[k] * x[k];

         if (r_best < r) r_best = r, delta = fabs(a[j]);

      }

      if (r_best < 0.001) r_best = 0.0;

      if (r_best == 0.0) goto done;

      xassert(delta > 0.0);

      /* try to increase violation by dividing delta by 2, 4, and 8,

         respectively */

      d_try[1] = delta / 2.0;

      d_try[2] = delta / 4.0;

      d_try[3] = delta / 8.0;

      for (j = 1; j <= 3; j++)

      {  /* construct c-MIR inequality */

         fail = cmir_ineq(n, a, b, u, cset, d_try[j], alpha, beta,

            gamma);

         if (fail) continue;

         /* compute violation */

         r = - (*beta) - (*gamma) * s;

         for (k = 1; k <= n; k++) r += alpha[k] * x[k];

         if (r_best < r) r_best = r, delta = d_try[j];

      }

      /* build subset of variables lying strictly between their bounds

         and order it by nondecreasing values of |x[j] - u[j]/2| */

      nv = 0;

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

      {  /* if x[j] is close to its bounds, skip it */

         eps = 1e-9 * (1.0 + fabs(u[j]));

         if (x[j] < eps || x[j] > u[j] - eps) continue;

         /* add x[j] to the subset */

         nv++;

         vset[nv].j = j;

         vset[nv].v = fabs(x[j] - 0.5 * u[j]);

      }

      qsort(&vset[1], nv, sizeof(struct vset), cmir_cmp);

      /* try to increase violation by successively complementing each

         variable in the subset */

      for (v = 1; v <= nv; v++)

      {  j = vset[v].j;

         /* replace x[j] by its complement or vice versa */

         cset[j] = (char)!cset[j];

         /* construct c-MIR inequality */

         fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);

         /* restore the variable */

         cset[j] = (char)!cset[j];

         /* do not replace the variable in case of failure */

         if (fail) continue;

         /* compute violation */

         r = - (*beta) - (*gamma) * s;

         for (k = 1; k <= n; k++) r += alpha[k] * x[k];

         if (r_best < r) r_best = r, cset[j] = (char)!cset[j];

      }

      /* construct the best c-MIR inequality chosen */

      fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);

      xassert(!fail);

done: /* free working arrays */

      xfree(cset);

      xfree(vset);

      /* return to the calling routine */

      return r_best;

}

static double generate(struct MIR *mir)

{     /* try to generate violated c-MIR cut for modified constraint */

      int m = mir->m;

      int n = mir->n;

      int j, k, kk, nint;

      double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma;

      ios_copy_vec(mir->cut_vec, mir->mod_vec);

      mir->cut_rhs = mir->mod_rhs;

      /* remove small terms, which can appear due to substitution of

         variable bounds */

      ios_clean_vec(mir->cut_vec, DBL_EPSILON);

#if _MIR_DEBUG

      ios_check_vec(mir->cut_vec);

#endif

      /* remove positive continuous terms to obtain MK relaxation */

      for (j = 1; j <= mir->cut_vec->nnz; j++)

      {  k = mir->cut_vec->ind[j];

         xassert(1 <= k && k <= m+n);

         if (!mir->isint[k] && mir->cut_vec->val[j] > 0.0)

            mir->cut_vec->val[j] = 0.0;

      }

      ios_clean_vec(mir->cut_vec, 0.0);

#if _MIR_DEBUG

      ios_check_vec(mir->cut_vec);

#endif

      /* move integer terms to the beginning of the sparse vector and

         determine the number of integer variables */

      nint = 0;

      for (j = 1; j <= mir->cut_vec->nnz; j++)

      {  k = mir->cut_vec->ind[j];

         xassert(1 <= k && k <= m+n);

         if (mir->isint[k])

         {  double temp;

            nint++;

            /* interchange elements [nint] and [j] */

            kk = mir->cut_vec->ind[nint];

            mir->cut_vec->pos[k] = nint;

            mir->cut_vec->pos[kk] = j;

            mir->cut_vec->ind[nint] = k;

            mir->cut_vec->ind[j] = kk;

            temp = mir->cut_vec->val[nint];

            mir->cut_vec->val[nint] = mir->cut_vec->val[j];

            mir->cut_vec->val[j] = temp;

         }

      }

#if _MIR_DEBUG

      ios_check_vec(mir->cut_vec);

#endif

      /* if there is no integer variable, nothing to generate */

      if (nint == 0) goto done;

      /* allocate working arrays */

      u = xcalloc(1+nint, sizeof(double));

      x = xcalloc(1+nint, sizeof(double));

      alpha = xcalloc(1+nint, sizeof(double));

      /* determine u and x */

      for (j = 1; j <= nint; j++)

      {  k = mir->cut_vec->ind[j];

         xassert(m+1 <= k && k <= m+n);

         xassert(mir->isint[k]);

         u[j] = mir->ub[k] - mir->lb[k];

         xassert(u[j] >= 1.0);

         if (mir->subst[k] == 'L')

            x[j] = mir->x[k] - mir->lb[k];

         else if (mir->subst[k] == 'U')

            x[j] = mir->ub[k] - mir->x[k];

         else

            xassert(k != k);

         xassert(x[j] >= -0.001);

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

      }

      /* compute s = - sum of continuous terms */

      s = 0.0;

      for (j = nint+1; j <= mir->cut_vec->nnz; j++)

      {  double x;

         k = mir->cut_vec->ind[j];

         xassert(1 <= k && k <= m+n);

         /* must be continuous */

         xassert(!mir->isint[k]);

         if (mir->subst[k] == 'L')

         {  xassert(mir->lb[k] != -DBL_MAX);

            kk = mir->vlb[k];

            if (kk == 0)

               x = mir->x[k] - mir->lb[k];

            else

               x = mir->x[k] - mir->lb[k] * mir->x[kk];

         }

         else if (mir->subst[k] == 'U')

         {  xassert(mir->ub[k] != +DBL_MAX);

            kk = mir->vub[k];

            if (kk == 0)

               x = mir->ub[k] - mir->x[k];

            else

               x = mir->ub[k] * mir->x[kk] - mir->x[k];

         }

         else

            xassert(k != k);

         xassert(x >= -0.001);

         if (x < 0.0) x = 0.0;

         s -= mir->cut_vec->val[j] * x;

      }

      xassert(s >= 0.0);

      /* apply heuristic to obtain most violated c-MIR inequality */

      b = mir->cut_rhs;

      r_best = cmir_sep(nint, mir->cut_vec->val, b, u, x, s, alpha,

         &beta, &gamma);

      if (r_best == 0.0) goto skip;

      xassert(r_best > 0.0);

      /* convert to raw cut */

      /* sum alpha[j] * x[j] <= beta + gamma * s */

      for (j = 1; j <= nint; j++)

         mir->cut_vec->val[j] = alpha[j];

      for (j = nint+1; j <= mir->cut_vec->nnz; j++)

      {  k = mir->cut_vec->ind[j];

         if (k <= m+n) mir->cut_vec->val[j] *= gamma;

      }

      mir->cut_rhs = beta;

#if _MIR_DEBUG

      ios_check_vec(mir->cut_vec);

#endif

skip: /* free working arrays */

      xfree(u);

      xfree(x);

      xfree(alpha);