alpar@1: /* glpios05.c (Gomory's mixed integer cut generator) */
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@1: *
alpar@1: *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@1: *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
alpar@1: *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@1: *  E-mail: <mao@gnu.org>.
alpar@1: *
alpar@1: *  GLPK is free software: you can redistribute it and/or modify it
alpar@1: *  under the terms of the GNU General Public License as published by
alpar@1: *  the Free Software Foundation, either version 3 of the License, or
alpar@1: *  (at your option) any later version.
alpar@1: *
alpar@1: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@1: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@1: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@1: *  License for more details.
alpar@1: *
alpar@1: *  You should have received a copy of the GNU General Public License
alpar@1: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@1: ***********************************************************************/
alpar@1: 
alpar@1: #include "glpios.h"
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  ios_gmi_gen - generate Gomory's mixed integer cuts.
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpios.h"
alpar@1: *  void ios_gmi_gen(glp_tree *tree, IOSPOOL *pool);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine ios_gmi_gen generates Gomory's mixed integer cuts for
alpar@1: *  the current point and adds them to the cut pool. */
alpar@1: 
alpar@1: #define MAXCUTS 50
alpar@1: /* maximal number of cuts to be generated for one round */
alpar@1: 
alpar@1: struct worka
alpar@1: {     /* Gomory's cut generator working area */
alpar@1:       int *ind; /* int ind[1+n]; */
alpar@1:       double *val; /* double val[1+n]; */
alpar@1:       double *phi; /* double phi[1+m+n]; */
alpar@1: };
alpar@1: 
alpar@1: #define f(x) ((x) - floor(x))
alpar@1: /* compute fractional part of x */
alpar@1: 
alpar@1: static void gen_cut(glp_tree *tree, struct worka *worka, int j)
alpar@1: {     /* this routine tries to generate Gomory's mixed integer cut for
alpar@1:          specified structural variable x[m+j] of integer kind, which is
alpar@1:          basic and has fractional value in optimal solution to current
alpar@1:          LP relaxation */
alpar@1:       glp_prob *mip = tree->mip;
alpar@1:       int m = mip->m;
alpar@1:       int n = mip->n;
alpar@1:       int *ind = worka->ind;
alpar@1:       double *val = worka->val;
alpar@1:       double *phi = worka->phi;
alpar@1:       int i, k, len, kind, stat;
alpar@1:       double lb, ub, alfa, beta, ksi, phi1, rhs;
alpar@1:       /* compute row of the simplex tableau, which (row) corresponds
alpar@1:          to specified basic variable xB[i] = x[m+j]; see (23) */
alpar@1:       len = glp_eval_tab_row(mip, m+j, ind, val);
alpar@1:       /* determine beta[i], which a value of xB[i] in optimal solution
alpar@1:          to current LP relaxation; note that this value is the same as
alpar@1:          if it would be computed with formula (27); it is assumed that
alpar@1:          beta[i] is fractional enough */
alpar@1:       beta = mip->col[j]->prim;
alpar@1:       /* compute cut coefficients phi and right-hand side rho, which
alpar@1:          correspond to formula (30); dense format is used, because rows
alpar@1:          of the simplex tableau is usually dense */
alpar@1:       for (k = 1; k <= m+n; k++) phi[k] = 0.0;
alpar@1:       rhs = f(beta); /* initial value of rho; see (28), (32) */
alpar@1:       for (j = 1; j <= len; j++)
alpar@1:       {  /* determine original number of non-basic variable xN[j] */
alpar@1:          k = ind[j];
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1:          /* determine the kind, bounds and current status of xN[j] in
alpar@1:             optimal solution to LP relaxation */
alpar@1:          if (k <= m)
alpar@1:          {  /* auxiliary variable */
alpar@1:             GLPROW *row = mip->row[k];
alpar@1:             kind = GLP_CV;
alpar@1:             lb = row->lb;
alpar@1:             ub = row->ub;
alpar@1:             stat = row->stat;
alpar@1:          }
alpar@1:          else
alpar@1:          {  /* structural variable */
alpar@1:             GLPCOL *col = mip->col[k-m];
alpar@1:             kind = col->kind;
alpar@1:             lb = col->lb;
alpar@1:             ub = col->ub;
alpar@1:             stat = col->stat;
alpar@1:          }
alpar@1:          /* xN[j] cannot be basic */
alpar@1:          xassert(stat != GLP_BS);
alpar@1:          /* determine row coefficient ksi[i,j] at xN[j]; see (23) */
alpar@1:          ksi = val[j];
alpar@1:          /* if ksi[i,j] is too large in the magnitude, do not generate
alpar@1:             the cut */
alpar@1:          if (fabs(ksi) > 1e+05) goto fini;
alpar@1:          /* if ksi[i,j] is too small in the magnitude, skip it */
alpar@1:          if (fabs(ksi) < 1e-10) goto skip;
alpar@1:          /* compute row coefficient alfa[i,j] at y[j]; see (26) */
alpar@1:          switch (stat)
alpar@1:          {  case GLP_NF:
alpar@1:                /* xN[j] is free (unbounded) having non-zero ksi[i,j];
alpar@1:                   do not generate the cut */
alpar@1:                goto fini;
alpar@1:             case GLP_NL:
alpar@1:                /* xN[j] has active lower bound */
alpar@1:                alfa = - ksi;
alpar@1:                break;
alpar@1:             case GLP_NU:
alpar@1:                /* xN[j] has active upper bound */
alpar@1:                alfa = + ksi;
alpar@1:                break;
alpar@1:             case GLP_NS:
alpar@1:                /* xN[j] is fixed; skip it */
alpar@1:                goto skip;
alpar@1:             default:
alpar@1:                xassert(stat != stat);
alpar@1:          }
alpar@1:          /* compute cut coefficient phi'[j] at y[j]; see (21), (28) */
alpar@1:          switch (kind)
alpar@1:          {  case GLP_IV:
alpar@1:                /* y[j] is integer */
alpar@1:                if (fabs(alfa - floor(alfa + 0.5)) < 1e-10)
alpar@1:                {  /* alfa[i,j] is close to nearest integer; skip it */
alpar@1:                   goto skip;
alpar@1:                }
alpar@1:                else if (f(alfa) <= f(beta))
alpar@1:                   phi1 = f(alfa);
alpar@1:                else
alpar@1:                   phi1 = (f(beta) / (1.0 - f(beta))) * (1.0 - f(alfa));
alpar@1:                break;
alpar@1:             case GLP_CV:
alpar@1:                /* y[j] is continuous */
alpar@1:                if (alfa >= 0.0)
alpar@1:                   phi1 = + alfa;
alpar@1:                else
alpar@1:                   phi1 = (f(beta) / (1.0 - f(beta))) * (- alfa);
alpar@1:                break;
alpar@1:             default:
alpar@1:                xassert(kind != kind);
alpar@1:          }
alpar@1:          /* compute cut coefficient phi[j] at xN[j] and update right-
alpar@1:             hand side rho; see (31), (32) */
alpar@1:          switch (stat)
alpar@1:          {  case GLP_NL:
alpar@1:                /* xN[j] has active lower bound */
alpar@1:                phi[k] = + phi1;
alpar@1:                rhs += phi1 * lb;
alpar@1:                break;
alpar@1:             case GLP_NU:
alpar@1:                /* xN[j] has active upper bound */
alpar@1:                phi[k] = - phi1;
alpar@1:                rhs -= phi1 * ub;
alpar@1:                break;
alpar@1:             default:
alpar@1:                xassert(stat != stat);
alpar@1:          }
alpar@1: skip:    ;
alpar@1:       }
alpar@1:       /* now the cut has the form sum_k phi[k] * x[k] >= rho, where cut
alpar@1:          coefficients are stored in the array phi in dense format;
alpar@1:          x[1,...,m] are auxiliary variables, x[m+1,...,m+n] are struc-
alpar@1:          tural variables; see (30) */
alpar@1:       /* eliminate auxiliary variables in order to express the cut only
alpar@1:          through structural variables; see (33) */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  GLPROW *row;
alpar@1:          GLPAIJ *aij;
alpar@1:          if (fabs(phi[i]) < 1e-10) continue;
alpar@1:          /* auxiliary variable x[i] has non-zero cut coefficient */
alpar@1:          row = mip->row[i];
alpar@1:          /* x[i] cannot be fixed */
alpar@1:          xassert(row->type != GLP_FX);
alpar@1:          /* substitute x[i] = sum_j a[i,j] * x[m+j] */
alpar@1:          for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@1:             phi[m+aij->col->j] += phi[i] * aij->val;
alpar@1:       }
alpar@1:       /* convert the final cut to sparse format and substitute fixed
alpar@1:          (structural) variables */
alpar@1:       len = 0;
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  GLPCOL *col;
alpar@1:          if (fabs(phi[m+j]) < 1e-10) continue;
alpar@1:          /* structural variable x[m+j] has non-zero cut coefficient */
alpar@1:          col = mip->col[j];
alpar@1:          if (col->type == GLP_FX)
alpar@1:          {  /* eliminate x[m+j] */
alpar@1:             rhs -= phi[m+j] * col->lb;
alpar@1:          }
alpar@1:          else
alpar@1:          {  len++;
alpar@1:             ind[len] = j;
alpar@1:             val[len] = phi[m+j];
alpar@1:          }
alpar@1:       }
alpar@1:       if (fabs(rhs) < 1e-12) rhs = 0.0;
alpar@1:       /* if the cut inequality seems to be badly scaled, reject it to
alpar@1:          avoid numeric difficulties */
alpar@1:       for (k = 1; k <= len; k++)
alpar@1:       {  if (fabs(val[k]) < 1e-03) goto fini;
alpar@1:          if (fabs(val[k]) > 1e+03) goto fini;
alpar@1:       }
alpar@1:       /* add the cut to the cut pool for further consideration */
alpar@1: #if 0
alpar@1:       ios_add_cut_row(tree, pool, GLP_RF_GMI, len, ind, val, GLP_LO,
alpar@1:          rhs);
alpar@1: #else
alpar@1:       glp_ios_add_row(tree, NULL, GLP_RF_GMI, 0, len, ind, val, GLP_LO,
alpar@1:          rhs);
alpar@1: #endif
alpar@1: fini: return;
alpar@1: }
alpar@1: 
alpar@1: struct var { int j; double f; };
alpar@1: 
alpar@1: static int fcmp(const void *p1, const void *p2)
alpar@1: {     const struct var *v1 = p1, *v2 = p2;
alpar@1:       if (v1->f > v2->f) return -1;
alpar@1:       if (v1->f < v2->f) return +1;
alpar@1:       return 0;
alpar@1: }
alpar@1: 
alpar@1: void ios_gmi_gen(glp_tree *tree)
alpar@1: {     /* main routine to generate Gomory's cuts */
alpar@1:       glp_prob *mip = tree->mip;
alpar@1:       int m = mip->m;
alpar@1:       int n = mip->n;
alpar@1:       struct var *var;
alpar@1:       int k, nv, j, size;
alpar@1:       struct worka _worka, *worka = &_worka;
alpar@1:       /* allocate working arrays */
alpar@1:       var = xcalloc(1+n, sizeof(struct var));
alpar@1:       worka->ind = xcalloc(1+n, sizeof(int));
alpar@1:       worka->val = xcalloc(1+n, sizeof(double));
alpar@1:       worka->phi = xcalloc(1+m+n, sizeof(double));
alpar@1:       /* build the list of integer structural variables, which are
alpar@1:          basic and have fractional value in optimal solution to current
alpar@1:          LP relaxation */
alpar@1:       nv = 0;
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  GLPCOL *col = mip->col[j];
alpar@1:          double frac;
alpar@1:          if (col->kind != GLP_IV) continue;
alpar@1:          if (col->type == GLP_FX) continue;
alpar@1:          if (col->stat != GLP_BS) continue;
alpar@1:          frac = f(col->prim);
alpar@1:          if (!(0.05 <= frac && frac <= 0.95)) continue;
alpar@1:          /* add variable to the list */
alpar@1:          nv++, var[nv].j = j, var[nv].f = frac;
alpar@1:       }
alpar@1:       /* order the list by descending fractionality */
alpar@1:       qsort(&var[1], nv, sizeof(struct var), fcmp);
alpar@1:       /* try to generate cuts by one for each variable in the list, but
alpar@1:          not more than MAXCUTS cuts */
alpar@1:       size = glp_ios_pool_size(tree);
alpar@1:       for (k = 1; k <= nv; k++)
alpar@1:       {  if (glp_ios_pool_size(tree) - size >= MAXCUTS) break;
alpar@1:          gen_cut(tree, worka, var[k].j);
alpar@1:       }
alpar@1:       /* free working arrays */
alpar@1:       xfree(var);
alpar@1:       xfree(worka->ind);
alpar@1:       xfree(worka->val);
alpar@1:       xfree(worka->phi);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /* eof */