alpar@1: /* glpios08.c (clique 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: static double get_row_lb(LPX *lp, int i)
alpar@1: {     /* this routine returns lower bound of row i or -DBL_MAX if the
alpar@1:          row has no lower bound */
alpar@1:       double lb;
alpar@1:       switch (lpx_get_row_type(lp, i))
alpar@1:       {  case LPX_FR:
alpar@1:          case LPX_UP:
alpar@1:             lb = -DBL_MAX;
alpar@1:             break;
alpar@1:          case LPX_LO:
alpar@1:          case LPX_DB:
alpar@1:          case LPX_FX:
alpar@1:             lb = lpx_get_row_lb(lp, i);
alpar@1:             break;
alpar@1:          default:
alpar@1:             xassert(lp != lp);
alpar@1:       }
alpar@1:       return lb;
alpar@1: }
alpar@1: 
alpar@1: static double get_row_ub(LPX *lp, int i)
alpar@1: {     /* this routine returns upper bound of row i or +DBL_MAX if the
alpar@1:          row has no upper bound */
alpar@1:       double ub;
alpar@1:       switch (lpx_get_row_type(lp, i))
alpar@1:       {  case LPX_FR:
alpar@1:          case LPX_LO:
alpar@1:             ub = +DBL_MAX;
alpar@1:             break;
alpar@1:          case LPX_UP:
alpar@1:          case LPX_DB:
alpar@1:          case LPX_FX:
alpar@1:             ub = lpx_get_row_ub(lp, i);
alpar@1:             break;
alpar@1:          default:
alpar@1:             xassert(lp != lp);
alpar@1:       }
alpar@1:       return ub;
alpar@1: }
alpar@1: 
alpar@1: static double get_col_lb(LPX *lp, int j)
alpar@1: {     /* this routine returns lower bound of column j or -DBL_MAX if
alpar@1:          the column has no lower bound */
alpar@1:       double lb;
alpar@1:       switch (lpx_get_col_type(lp, j))
alpar@1:       {  case LPX_FR:
alpar@1:          case LPX_UP:
alpar@1:             lb = -DBL_MAX;
alpar@1:             break;
alpar@1:          case LPX_LO:
alpar@1:          case LPX_DB:
alpar@1:          case LPX_FX:
alpar@1:             lb = lpx_get_col_lb(lp, j);
alpar@1:             break;
alpar@1:          default:
alpar@1:             xassert(lp != lp);
alpar@1:       }
alpar@1:       return lb;
alpar@1: }
alpar@1: 
alpar@1: static double get_col_ub(LPX *lp, int j)
alpar@1: {     /* this routine returns upper bound of column j or +DBL_MAX if
alpar@1:          the column has no upper bound */
alpar@1:       double ub;
alpar@1:       switch (lpx_get_col_type(lp, j))
alpar@1:       {  case LPX_FR:
alpar@1:          case LPX_LO:
alpar@1:             ub = +DBL_MAX;
alpar@1:             break;
alpar@1:          case LPX_UP:
alpar@1:          case LPX_DB:
alpar@1:          case LPX_FX:
alpar@1:             ub = lpx_get_col_ub(lp, j);
alpar@1:             break;
alpar@1:          default:
alpar@1:             xassert(lp != lp);
alpar@1:       }
alpar@1:       return ub;
alpar@1: }
alpar@1: 
alpar@1: static int is_binary(LPX *lp, int j)
alpar@1: {     /* this routine checks if variable x[j] is binary */
alpar@1:       return
alpar@1:          lpx_get_col_kind(lp, j) == LPX_IV &&
alpar@1:          lpx_get_col_type(lp, j) == LPX_DB &&
alpar@1:          lpx_get_col_lb(lp, j) == 0.0 && lpx_get_col_ub(lp, j) == 1.0;
alpar@1: }
alpar@1: 
alpar@1: static double eval_lf_min(LPX *lp, int len, int ind[], double val[])
alpar@1: {     /* this routine computes the minimum of a specified linear form
alpar@1: 
alpar@1:             sum a[j]*x[j]
alpar@1:              j
alpar@1: 
alpar@1:          using the formula:
alpar@1: 
alpar@1:             min =   sum   a[j]*lb[j] +   sum   a[j]*ub[j],
alpar@1:                   j in J+              j in J-
alpar@1: 
alpar@1:          where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j]
alpar@1:          are lower and upper bound of variable x[j], resp. */
alpar@1:       int j, t;
alpar@1:       double lb, ub, sum;
alpar@1:       sum = 0.0;
alpar@1:       for (t = 1; t <= len; t++)
alpar@1:       {  j = ind[t];
alpar@1:          if (val[t] > 0.0)
alpar@1:          {  lb = get_col_lb(lp, j);
alpar@1:             if (lb == -DBL_MAX)
alpar@1:             {  sum = -DBL_MAX;
alpar@1:                break;
alpar@1:             }
alpar@1:             sum += val[t] * lb;
alpar@1:          }
alpar@1:          else if (val[t] < 0.0)
alpar@1:          {  ub = get_col_ub(lp, j);
alpar@1:             if (ub == +DBL_MAX)
alpar@1:             {  sum = -DBL_MAX;
alpar@1:                break;
alpar@1:             }
alpar@1:             sum += val[t] * ub;
alpar@1:          }
alpar@1:          else
alpar@1:             xassert(val != val);
alpar@1:       }
alpar@1:       return sum;
alpar@1: }
alpar@1: 
alpar@1: static double eval_lf_max(LPX *lp, int len, int ind[], double val[])
alpar@1: {     /* this routine computes the maximum of a specified linear form
alpar@1: 
alpar@1:             sum a[j]*x[j]
alpar@1:              j
alpar@1: 
alpar@1:          using the formula:
alpar@1: 
alpar@1:             max =   sum   a[j]*ub[j] +   sum   a[j]*lb[j],
alpar@1:                   j in J+              j in J-
alpar@1: 
alpar@1:          where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j]
alpar@1:          are lower and upper bound of variable x[j], resp. */
alpar@1:       int j, t;
alpar@1:       double lb, ub, sum;
alpar@1:       sum = 0.0;
alpar@1:       for (t = 1; t <= len; t++)
alpar@1:       {  j = ind[t];
alpar@1:          if (val[t] > 0.0)
alpar@1:          {  ub = get_col_ub(lp, j);
alpar@1:             if (ub == +DBL_MAX)
alpar@1:             {  sum = +DBL_MAX;
alpar@1:                break;
alpar@1:             }
alpar@1:             sum += val[t] * ub;
alpar@1:          }
alpar@1:          else if (val[t] < 0.0)
alpar@1:          {  lb = get_col_lb(lp, j);
alpar@1:             if (lb == -DBL_MAX)
alpar@1:             {  sum = +DBL_MAX;
alpar@1:                break;
alpar@1:             }
alpar@1:             sum += val[t] * lb;
alpar@1:          }
alpar@1:          else
alpar@1:             xassert(val != val);
alpar@1:       }
alpar@1:       return sum;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: -- probing - determine logical relation between binary variables.
alpar@1: --
alpar@1: -- This routine tentatively sets a binary variable to 0 and then to 1
alpar@1: -- and examines whether another binary variable is caused to be fixed.
alpar@1: --
alpar@1: -- The examination is based only on one row (constraint), which is the
alpar@1: -- following:
alpar@1: --
alpar@1: --    L <= sum a[j]*x[j] <= U.                                       (1)
alpar@1: --          j
alpar@1: --
alpar@1: -- Let x[p] be a probing variable, x[q] be an examined variable. Then
alpar@1: -- (1) can be written as:
alpar@1: --
alpar@1: --    L <=   sum  a[j]*x[j] + a[p]*x[p] + a[q]*x[q] <= U,            (2)
alpar@1: --         j in J'
alpar@1: --
alpar@1: -- where J' = {j: j != p and j != q}.
alpar@1: --
alpar@1: -- Let
alpar@1: --
alpar@1: --    L' = L - a[p]*x[p],                                            (3)
alpar@1: --
alpar@1: --    U' = U - a[p]*x[p],                                            (4)
alpar@1: --
alpar@1: -- where x[p] is assumed to be fixed at 0 or 1. So (2) can be rewritten
alpar@1: -- as follows:
alpar@1: --
alpar@1: --    L' <=   sum  a[j]*x[j] + a[q]*x[q] <= U',                      (5)
alpar@1: --          j in J'
alpar@1: --
alpar@1: -- from where we have:
alpar@1: --
alpar@1: --    L' -  sum  a[j]*x[j] <= a[q]*x[q] <= U' -  sum  a[j]*x[j].     (6)
alpar@1: --        j in J'                              j in J'
alpar@1: --
alpar@1: -- Thus,
alpar@1: --
alpar@1: --    min a[q]*x[q] = L' - MAX,                                      (7)
alpar@1: --
alpar@1: --    max a[q]*x[q] = U' - MIN,                                      (8)
alpar@1: --
alpar@1: -- where
alpar@1: --
alpar@1: --    MIN = min  sum  a[j]*x[j],                                     (9)
alpar@1: --             j in J'
alpar@1: --
alpar@1: --    MAX = max  sum  a[j]*x[j].                                    (10)
alpar@1: --             j in J'
alpar@1: --
alpar@1: -- Formulae (7) and (8) allows determining implied lower and upper
alpar@1: -- bounds of x[q].
alpar@1: --
alpar@1: -- Parameters len, val, L and U specify the constraint (1).
alpar@1: --
alpar@1: -- Parameters lf_min and lf_max specify implied lower and upper bounds
alpar@1: -- of the linear form (1). It is assumed that these bounds are computed
alpar@1: -- with the routines eval_lf_min and eval_lf_max (see above).
alpar@1: --
alpar@1: -- Parameter p specifies the probing variable x[p], which is set to 0
alpar@1: -- (if set is 0) or to 1 (if set is 1).
alpar@1: --
alpar@1: -- Parameter q specifies the examined variable x[q].
alpar@1: --
alpar@1: -- On exit the routine returns one of the following codes:
alpar@1: --
alpar@1: -- 0 - there is no logical relation between x[p] and x[q];
alpar@1: -- 1 - x[q] can take only on value 0;
alpar@1: -- 2 - x[q] can take only on value 1. */
alpar@1: 
alpar@1: static int probing(int len, double val[], double L, double U,
alpar@1:       double lf_min, double lf_max, int p, int set, int q)
alpar@1: {     double temp;
alpar@1:       xassert(1 <= p && p < q && q <= len);
alpar@1:       /* compute L' (3) */
alpar@1:       if (L != -DBL_MAX && set) L -= val[p];
alpar@1:       /* compute U' (4) */
alpar@1:       if (U != +DBL_MAX && set) U -= val[p];
alpar@1:       /* compute MIN (9) */
alpar@1:       if (lf_min != -DBL_MAX)
alpar@1:       {  if (val[p] < 0.0) lf_min -= val[p];
alpar@1:          if (val[q] < 0.0) lf_min -= val[q];
alpar@1:       }
alpar@1:       /* compute MAX (10) */
alpar@1:       if (lf_max != +DBL_MAX)
alpar@1:       {  if (val[p] > 0.0) lf_max -= val[p];
alpar@1:          if (val[q] > 0.0) lf_max -= val[q];
alpar@1:       }
alpar@1:       /* compute implied lower bound of x[q]; see (7), (8) */
alpar@1:       if (val[q] > 0.0)
alpar@1:       {  if (L == -DBL_MAX || lf_max == +DBL_MAX)
alpar@1:             temp = -DBL_MAX;
alpar@1:          else
alpar@1:             temp = (L - lf_max) / val[q];
alpar@1:       }
alpar@1:       else
alpar@1:       {  if (U == +DBL_MAX || lf_min == -DBL_MAX)
alpar@1:             temp = -DBL_MAX;
alpar@1:          else
alpar@1:             temp = (U - lf_min) / val[q];
alpar@1:       }
alpar@1:       if (temp > 0.001) return 2;
alpar@1:       /* compute implied upper bound of x[q]; see (7), (8) */
alpar@1:       if (val[q] > 0.0)
alpar@1:       {  if (U == +DBL_MAX || lf_min == -DBL_MAX)
alpar@1:             temp = +DBL_MAX;
alpar@1:          else
alpar@1:             temp = (U - lf_min) / val[q];
alpar@1:       }
alpar@1:       else
alpar@1:       {  if (L == -DBL_MAX || lf_max == +DBL_MAX)
alpar@1:             temp = +DBL_MAX;
alpar@1:          else
alpar@1:             temp = (L - lf_max) / val[q];
alpar@1:       }
alpar@1:       if (temp < 0.999) return 1;
alpar@1:       /* there is no logical relation between x[p] and x[q] */
alpar@1:       return 0;
alpar@1: }
alpar@1: 
alpar@1: struct COG
alpar@1: {     /* conflict graph; it represents logical relations between binary
alpar@1:          variables and has a vertex for each binary variable and its
alpar@1:          complement, and an edge between two vertices when at most one
alpar@1:          of the variables represented by the vertices can equal one in
alpar@1:          an optimal solution */
alpar@1:       int n;
alpar@1:       /* number of variables */
alpar@1:       int nb;
alpar@1:       /* number of binary variables represented in the graph (note that
alpar@1:          not all binary variables can be represented); vertices which
alpar@1:          correspond to binary variables have numbers 1, ..., nb while
alpar@1:          vertices which correspond to complements of binary variables
alpar@1:          have numbers nb+1, ..., nb+nb */
alpar@1:       int ne;
alpar@1:       /* number of edges in the graph */
alpar@1:       int *vert; /* int vert[1+n]; */
alpar@1:       /* if x[j] is a binary variable represented in the graph, vert[j]
alpar@1:          is the vertex number corresponding to x[j]; otherwise vert[j]
alpar@1:          is zero */
alpar@1:       int *orig; /* int list[1:nb]; */
alpar@1:       /* if vert[j] = k > 0, then orig[k] = j */
alpar@1:       unsigned char *a;
alpar@1:       /* adjacency matrix of the graph having 2*nb rows and columns;
alpar@1:          only strict lower triangle is stored in dense packed form */
alpar@1: };
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: -- lpx_create_cog - create the conflict graph.
alpar@1: --
alpar@1: -- SYNOPSIS
alpar@1: --
alpar@1: -- #include "glplpx.h"
alpar@1: -- void *lpx_create_cog(LPX *lp);
alpar@1: --
alpar@1: -- DESCRIPTION
alpar@1: --
alpar@1: -- The routine lpx_create_cog creates the conflict graph for a given
alpar@1: -- problem instance.
alpar@1: --
alpar@1: -- RETURNS
alpar@1: --
alpar@1: -- If the graph has been created, the routine returns a pointer to it.
alpar@1: -- Otherwise the routine returns NULL. */
alpar@1: 
alpar@1: #define MAX_NB 4000
alpar@1: #define MAX_ROW_LEN 500
alpar@1: 
alpar@1: static void lpx_add_cog_edge(void *_cog, int i, int j);
alpar@1: 
alpar@1: static void *lpx_create_cog(LPX *lp)
alpar@1: {     struct COG *cog = NULL;
alpar@1:       int m, n, nb, i, j, p, q, len, *ind, *vert, *orig;
alpar@1:       double L, U, lf_min, lf_max, *val;
alpar@1:       xprintf("Creating the conflict graph...\n");
alpar@1:       m = lpx_get_num_rows(lp);
alpar@1:       n = lpx_get_num_cols(lp);
alpar@1:       /* determine which binary variables should be included in the
alpar@1:          conflict graph */
alpar@1:       nb = 0;
alpar@1:       vert = xcalloc(1+n, sizeof(int));
alpar@1:       for (j = 1; j <= n; j++) vert[j] = 0;
alpar@1:       orig = xcalloc(1+n, sizeof(int));
alpar@1:       ind = xcalloc(1+n, sizeof(int));
alpar@1:       val = xcalloc(1+n, sizeof(double));
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  L = get_row_lb(lp, i);
alpar@1:          U = get_row_ub(lp, i);
alpar@1:          if (L == -DBL_MAX && U == +DBL_MAX) continue;
alpar@1:          len = lpx_get_mat_row(lp, i, ind, val);
alpar@1:          if (len > MAX_ROW_LEN) continue;
alpar@1:          lf_min = eval_lf_min(lp, len, ind, val);
alpar@1:          lf_max = eval_lf_max(lp, len, ind, val);
alpar@1:          for (p = 1; p <= len; p++)
alpar@1:          {  if (!is_binary(lp, ind[p])) continue;
alpar@1:             for (q = p+1; q <= len; q++)
alpar@1:             {  if (!is_binary(lp, ind[q])) continue;
alpar@1:                if (probing(len, val, L, U, lf_min, lf_max, p, 0, q) ||
alpar@1:                    probing(len, val, L, U, lf_min, lf_max, p, 1, q))
alpar@1:                {  /* there is a logical relation */
alpar@1:                   /* include the first variable in the graph */
alpar@1:                   j = ind[p];
alpar@1:                   if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j;
alpar@1:                   /* incude the second variable in the graph */
alpar@1:                   j = ind[q];
alpar@1:                   if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j;
alpar@1:                }
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       /* if the graph is either empty or has too many vertices, do not
alpar@1:          create it */
alpar@1:       if (nb == 0 || nb > MAX_NB)
alpar@1:       {  xprintf("The conflict graph is either empty or too big\n");
alpar@1:          xfree(vert);
alpar@1:          xfree(orig);
alpar@1:          goto done;
alpar@1:       }
alpar@1:       /* create the conflict graph */
alpar@1:       cog = xmalloc(sizeof(struct COG));
alpar@1:       cog->n = n;
alpar@1:       cog->nb = nb;
alpar@1:       cog->ne = 0;
alpar@1:       cog->vert = vert;
alpar@1:       cog->orig = orig;
alpar@1:       len = nb + nb; /* number of vertices */
alpar@1:       len = (len * (len - 1)) / 2; /* number of entries in triangle */
alpar@1:       len = (len + (CHAR_BIT - 1)) / CHAR_BIT; /* bytes needed */
alpar@1:       cog->a = xmalloc(len);
alpar@1:       memset(cog->a, 0, len);
alpar@1:       for (j = 1; j <= nb; j++)
alpar@1:       {  /* add edge between variable and its complement */
alpar@1:          lpx_add_cog_edge(cog, +orig[j], -orig[j]);
alpar@1:       }
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  L = get_row_lb(lp, i);
alpar@1:          U = get_row_ub(lp, i);
alpar@1:          if (L == -DBL_MAX && U == +DBL_MAX) continue;
alpar@1:          len = lpx_get_mat_row(lp, i, ind, val);
alpar@1:          if (len > MAX_ROW_LEN) continue;
alpar@1:          lf_min = eval_lf_min(lp, len, ind, val);
alpar@1:          lf_max = eval_lf_max(lp, len, ind, val);
alpar@1:          for (p = 1; p <= len; p++)
alpar@1:          {  if (!is_binary(lp, ind[p])) continue;
alpar@1:             for (q = p+1; q <= len; q++)
alpar@1:             {  if (!is_binary(lp, ind[q])) continue;
alpar@1:                /* set x[p] to 0 and examine x[q] */
alpar@1:                switch (probing(len, val, L, U, lf_min, lf_max, p, 0, q))
alpar@1:                {  case 0:
alpar@1:                      /* no logical relation */
alpar@1:                      break;
alpar@1:                   case 1:
alpar@1:                      /* x[p] = 0 implies x[q] = 0 */
alpar@1:                      lpx_add_cog_edge(cog, -ind[p], +ind[q]);
alpar@1:                      break;
alpar@1:                   case 2:
alpar@1:                      /* x[p] = 0 implies x[q] = 1 */
alpar@1:                      lpx_add_cog_edge(cog, -ind[p], -ind[q]);
alpar@1:                      break;
alpar@1:                   default:
alpar@1:                      xassert(lp != lp);
alpar@1:                }
alpar@1:                /* set x[p] to 1 and examine x[q] */
alpar@1:                switch (probing(len, val, L, U, lf_min, lf_max, p, 1, q))
alpar@1:                {  case 0:
alpar@1:                      /* no logical relation */
alpar@1:                      break;
alpar@1:                   case 1:
alpar@1:                      /* x[p] = 1 implies x[q] = 0 */
alpar@1:                      lpx_add_cog_edge(cog, +ind[p], +ind[q]);
alpar@1:                      break;
alpar@1:                   case 2:
alpar@1:                      /* x[p] = 1 implies x[q] = 1 */
alpar@1:                      lpx_add_cog_edge(cog, +ind[p], -ind[q]);
alpar@1:                      break;
alpar@1:                   default:
alpar@1:                      xassert(lp != lp);
alpar@1:                }
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       xprintf("The conflict graph has 2*%d vertices and %d edges\n",
alpar@1:          cog->nb, cog->ne);
alpar@1: done: xfree(ind);
alpar@1:       xfree(val);
alpar@1:       return cog;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: -- lpx_add_cog_edge - add edge to the conflict graph.
alpar@1: --
alpar@1: -- SYNOPSIS
alpar@1: --
alpar@1: -- #include "glplpx.h"
alpar@1: -- void lpx_add_cog_edge(void *cog, int i, int j);
alpar@1: --
alpar@1: -- DESCRIPTION
alpar@1: --
alpar@1: -- The routine lpx_add_cog_edge adds an edge to the conflict graph.
alpar@1: -- The edge connects x[i] (if i > 0) or its complement (if i < 0) and
alpar@1: -- x[j] (if j > 0) or its complement (if j < 0), where i and j are
alpar@1: -- original ordinal numbers of corresponding variables. */
alpar@1: 
alpar@1: static void lpx_add_cog_edge(void *_cog, int i, int j)
alpar@1: {     struct COG *cog = _cog;
alpar@1:       int k;
alpar@1:       xassert(i != j);
alpar@1:       /* determine indices of corresponding vertices */
alpar@1:       if (i > 0)
alpar@1:       {  xassert(1 <= i && i <= cog->n);
alpar@1:          i = cog->vert[i];
alpar@1:          xassert(i != 0);
alpar@1:       }
alpar@1:       else
alpar@1:       {  i = -i;
alpar@1:          xassert(1 <= i && i <= cog->n);
alpar@1:          i = cog->vert[i];
alpar@1:          xassert(i != 0);
alpar@1:          i += cog->nb;
alpar@1:       }
alpar@1:       if (j > 0)
alpar@1:       {  xassert(1 <= j && j <= cog->n);
alpar@1:          j = cog->vert[j];
alpar@1:          xassert(j != 0);
alpar@1:       }
alpar@1:       else
alpar@1:       {  j = -j;
alpar@1:          xassert(1 <= j && j <= cog->n);
alpar@1:          j = cog->vert[j];
alpar@1:          xassert(j != 0);
alpar@1:          j += cog->nb;
alpar@1:       }
alpar@1:       /* only lower triangle is stored, so we need i > j */
alpar@1:       if (i < j) k = i, i = j, j = k;
alpar@1:       k = ((i - 1) * (i - 2)) / 2 + (j - 1);
alpar@1:       cog->a[k / CHAR_BIT] |=
alpar@1:          (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
alpar@1:       cog->ne++;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: -- MAXIMUM WEIGHT CLIQUE
alpar@1: --
alpar@1: -- Two subroutines sub() and wclique() below are intended to find a
alpar@1: -- maximum weight clique in a given undirected graph. These subroutines
alpar@1: -- are slightly modified version of the program WCLIQUE developed by
alpar@1: -- Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based
alpar@1: -- on ideas from the article "P. R. J. Ostergard, A new algorithm for
alpar@1: -- the maximum-weight clique problem, submitted for publication", which
alpar@1: -- in turn is a generalization of the algorithm for unweighted graphs
alpar@1: -- presented in "P. R. J. Ostergard, A fast algorithm for the maximum
alpar@1: -- clique problem, submitted for publication".
alpar@1: --
alpar@1: -- USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE. */
alpar@1: 
alpar@1: struct dsa
alpar@1: {     /* dynamic storage area */
alpar@1:       int n;
alpar@1:       /* number of vertices */
alpar@1:       int *wt; /* int wt[0:n-1]; */
alpar@1:       /* weights */
alpar@1:       unsigned char *a;
alpar@1:       /* adjacency matrix (packed lower triangle without main diag.) */
alpar@1:       int record;
alpar@1:       /* weight of best clique */
alpar@1:       int rec_level;
alpar@1:       /* number of vertices in best clique */
alpar@1:       int *rec; /* int rec[0:n-1]; */
alpar@1:       /* best clique so far */
alpar@1:       int *clique; /* int clique[0:n-1]; */
alpar@1:       /* table for pruning */
alpar@1:       int *set; /* int set[0:n-1]; */
alpar@1:       /* current clique */
alpar@1: };
alpar@1: 
alpar@1: #define n         (dsa->n)
alpar@1: #define wt        (dsa->wt)
alpar@1: #define a         (dsa->a)
alpar@1: #define record    (dsa->record)
alpar@1: #define rec_level (dsa->rec_level)
alpar@1: #define rec       (dsa->rec)
alpar@1: #define clique    (dsa->clique)
alpar@1: #define set       (dsa->set)
alpar@1: 
alpar@1: #if 0
alpar@1: static int is_edge(struct dsa *dsa, int i, int j)
alpar@1: {     /* if there is arc (i,j), the routine returns true; otherwise
alpar@1:          false; 0 <= i, j < n */
alpar@1:       int k;
alpar@1:       xassert(0 <= i && i < n);
alpar@1:       xassert(0 <= j && j < n);
alpar@1:       if (i == j) return 0;
alpar@1:       if (i < j) k = i, i = j, j = k;
alpar@1:       k = (i * (i - 1)) / 2 + j;
alpar@1:       return a[k / CHAR_BIT] &
alpar@1:          (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
alpar@1: }
alpar@1: #else
alpar@1: #define is_edge(dsa, i, j) ((i) == (j) ? 0 : \
alpar@1:       (i) > (j) ? is_edge1(i, j) : is_edge1(j, i))
alpar@1: #define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j))
alpar@1: #define is_edge2(k) (a[(k) / CHAR_BIT] & \
alpar@1:       (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT)))
alpar@1: #endif
alpar@1: 
alpar@1: static void sub(struct dsa *dsa, int ct, int table[], int level,
alpar@1:       int weight, int l_weight)
alpar@1: {     int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable;
alpar@1:       newtable = xcalloc(n, sizeof(int));
alpar@1:       if (ct <= 0)
alpar@1:       {  /* 0 or 1 elements left; include these */
alpar@1:          if (ct == 0)
alpar@1:          {  set[level++] = table[0];
alpar@1:             weight += l_weight;
alpar@1:          }
alpar@1:          if (weight > record)
alpar@1:          {  record = weight;
alpar@1:             rec_level = level;
alpar@1:             for (i = 0; i < level; i++) rec[i] = set[i];
alpar@1:          }
alpar@1:          goto done;
alpar@1:       }
alpar@1:       for (i = ct; i >= 0; i--)
alpar@1:       {  if ((level == 0) && (i < ct)) goto done;
alpar@1:          k = table[i];
alpar@1:          if ((level > 0) && (clique[k] <= (record - weight)))
alpar@1:             goto done; /* prune */
alpar@1:          set[level] = k;
alpar@1:          curr_weight = weight + wt[k];
alpar@1:          l_weight -= wt[k];
alpar@1:          if (l_weight <= (record - curr_weight))
alpar@1:             goto done; /* prune */
alpar@1:          p1 = newtable;
alpar@1:          p2 = table;
alpar@1:          left_weight = 0;
alpar@1:          while (p2 < table + i)
alpar@1:          {  j = *p2++;
alpar@1:             if (is_edge(dsa, j, k))
alpar@1:             {  *p1++ = j;
alpar@1:                left_weight += wt[j];
alpar@1:             }
alpar@1:          }
alpar@1:          if (left_weight <= (record - curr_weight)) continue;
alpar@1:          sub(dsa, p1 - newtable - 1, newtable, level + 1, curr_weight,
alpar@1:             left_weight);
alpar@1:       }
alpar@1: done: xfree(newtable);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: static int wclique(int _n, int w[], unsigned char _a[], int sol[])
alpar@1: {     struct dsa _dsa, *dsa = &_dsa;
alpar@1:       int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos;
alpar@1:       glp_long timer;
alpar@1:       n = _n;
alpar@1:       wt = &w[1];
alpar@1:       a = _a;
alpar@1:       record = 0;
alpar@1:       rec_level = 0;
alpar@1:       rec = &sol[1];
alpar@1:       clique = xcalloc(n, sizeof(int));
alpar@1:       set = xcalloc(n, sizeof(int));
alpar@1:       used = xcalloc(n, sizeof(int));
alpar@1:       nwt = xcalloc(n, sizeof(int));
alpar@1:       pos = xcalloc(n, sizeof(int));
alpar@1:       /* start timer */
alpar@1:       timer = xtime();
alpar@1:       /* order vertices */
alpar@1:       for (i = 0; i < n; i++)
alpar@1:       {  nwt[i] = 0;
alpar@1:          for (j = 0; j < n; j++)
alpar@1:             if (is_edge(dsa, i, j)) nwt[i] += wt[j];
alpar@1:       }
alpar@1:       for (i = 0; i < n; i++)
alpar@1:          used[i] = 0;
alpar@1:       for (i = n-1; i >= 0; i--)
alpar@1:       {  max_wt = -1;
alpar@1:          max_nwt = -1;
alpar@1:          for (j = 0; j < n; j++)
alpar@1:          {  if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt
alpar@1:                && nwt[j] > max_nwt)))
alpar@1:             {  max_wt = wt[j];
alpar@1:                max_nwt = nwt[j];
alpar@1:                p = j;
alpar@1:             }
alpar@1:          }
alpar@1:          pos[i] = p;
alpar@1:          used[p] = 1;
alpar@1:          for (j = 0; j < n; j++)
alpar@1:             if ((!used[j]) && (j != p) && (is_edge(dsa, p, j)))
alpar@1:                nwt[j] -= wt[p];
alpar@1:       }
alpar@1:       /* main routine */
alpar@1:       wth = 0;
alpar@1:       for (i = 0; i < n; i++)
alpar@1:       {  wth += wt[pos[i]];
alpar@1:          sub(dsa, i, pos, 0, 0, wth);
alpar@1:          clique[pos[i]] = record;
alpar@1: #if 0
alpar@1:          if (utime() >= timer + 5.0)
alpar@1: #else
alpar@1:          if (xdifftime(xtime(), timer) >= 5.0 - 0.001)
alpar@1: #endif
alpar@1:          {  /* print current record and reset timer */
alpar@1:             xprintf("level = %d (%d); best = %d\n", i+1, n, record);
alpar@1: #if 0
alpar@1:             timer = utime();
alpar@1: #else
alpar@1:             timer = xtime();
alpar@1: #endif
alpar@1:          }
alpar@1:       }
alpar@1:       xfree(clique);
alpar@1:       xfree(set);
alpar@1:       xfree(used);
alpar@1:       xfree(nwt);
alpar@1:       xfree(pos);
alpar@1:       /* return the solution found */
alpar@1:       for (i = 1; i <= rec_level; i++) sol[i]++;
alpar@1:       return rec_level;
alpar@1: }
alpar@1: 
alpar@1: #undef n
alpar@1: #undef wt
alpar@1: #undef a
alpar@1: #undef record
alpar@1: #undef rec_level
alpar@1: #undef rec
alpar@1: #undef clique
alpar@1: #undef set
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: -- lpx_clique_cut - generate cluque cut.
alpar@1: --
alpar@1: -- SYNOPSIS
alpar@1: --
alpar@1: -- #include "glplpx.h"
alpar@1: -- int lpx_clique_cut(LPX *lp, void *cog, int ind[], double val[]);
alpar@1: --
alpar@1: -- DESCRIPTION
alpar@1: --
alpar@1: -- The routine lpx_clique_cut generates a clique cut using the conflict
alpar@1: -- graph specified by the parameter cog.
alpar@1: --
alpar@1: -- If a violated clique cut has been found, it has the following form:
alpar@1: --
alpar@1: --    sum{j in J} a[j]*x[j] <= b.
alpar@1: --
alpar@1: -- Variable indices j in J are stored in elements ind[1], ..., ind[len]
alpar@1: -- while corresponding constraint coefficients are stored in elements
alpar@1: -- val[1], ..., val[len], where len is returned on exit. The right-hand
alpar@1: -- side b is stored in element val[0].
alpar@1: --
alpar@1: -- RETURNS
alpar@1: --
alpar@1: -- If the cutting plane has been successfully generated, the routine
alpar@1: -- returns 1 <= len <= n, which is the number of non-zero coefficients
alpar@1: -- in the inequality constraint. Otherwise, the routine returns zero. */
alpar@1: 
alpar@1: static int lpx_clique_cut(LPX *lp, void *_cog, int ind[], double val[])
alpar@1: {     struct COG *cog = _cog;
alpar@1:       int n = lpx_get_num_cols(lp);
alpar@1:       int j, t, v, card, temp, len = 0, *w, *sol;
alpar@1:       double x, sum, b, *vec;
alpar@1:       /* allocate working arrays */
alpar@1:       w = xcalloc(1 + 2 * cog->nb, sizeof(int));
alpar@1:       sol = xcalloc(1 + 2 * cog->nb, sizeof(int));
alpar@1:       vec = xcalloc(1+n, sizeof(double));
alpar@1:       /* assign weights to vertices of the conflict graph */
alpar@1:       for (t = 1; t <= cog->nb; t++)
alpar@1:       {  j = cog->orig[t];
alpar@1:          x = lpx_get_col_prim(lp, j);
alpar@1:          temp = (int)(100.0 * x + 0.5);
alpar@1:          if (temp < 0) temp = 0;
alpar@1:          if (temp > 100) temp = 100;
alpar@1:          w[t] = temp;
alpar@1:          w[cog->nb + t] = 100 - temp;
alpar@1:       }
alpar@1:       /* find a clique of maximum weight */
alpar@1:       card = wclique(2 * cog->nb, w, cog->a, sol);
alpar@1:       /* compute the clique weight for unscaled values */
alpar@1:       sum = 0.0;
alpar@1:       for ( t = 1; t <= card; t++)
alpar@1:       {  v = sol[t];
alpar@1:          xassert(1 <= v && v <= 2 * cog->nb);
alpar@1:          if (v <= cog->nb)
alpar@1:          {  /* vertex v corresponds to binary variable x[j] */
alpar@1:             j = cog->orig[v];
alpar@1:             x = lpx_get_col_prim(lp, j);
alpar@1:             sum += x;
alpar@1:          }
alpar@1:          else
alpar@1:          {  /* vertex v corresponds to the complement of x[j] */
alpar@1:             j = cog->orig[v - cog->nb];
alpar@1:             x = lpx_get_col_prim(lp, j);
alpar@1:             sum += 1.0 - x;
alpar@1:          }
alpar@1:       }
alpar@1:       /* if the sum of binary variables and their complements in the
alpar@1:          clique greater than 1, the clique cut is violated */
alpar@1:       if (sum >= 1.01)
alpar@1:       {  /* construct the inquality */
alpar@1:          for (j = 1; j <= n; j++) vec[j] = 0;
alpar@1:          b = 1.0;
alpar@1:          for (t = 1; t <= card; t++)
alpar@1:          {  v = sol[t];
alpar@1:             if (v <= cog->nb)
alpar@1:             {  /* vertex v corresponds to binary variable x[j] */
alpar@1:                j = cog->orig[v];
alpar@1:                xassert(1 <= j && j <= n);
alpar@1:                vec[j] += 1.0;
alpar@1:             }
alpar@1:             else
alpar@1:             {  /* vertex v corresponds to the complement of x[j] */
alpar@1:                j = cog->orig[v - cog->nb];
alpar@1:                xassert(1 <= j && j <= n);
alpar@1:                vec[j] -= 1.0;
alpar@1:                b -= 1.0;
alpar@1:             }
alpar@1:          }
alpar@1:          xassert(len == 0);
alpar@1:          for (j = 1; j <= n; j++)
alpar@1:          {  if (vec[j] != 0.0)
alpar@1:             {  len++;
alpar@1:                ind[len] = j, val[len] = vec[j];
alpar@1:             }
alpar@1:          }
alpar@1:          ind[0] = 0, val[0] = b;
alpar@1:       }
alpar@1:       /* free working arrays */
alpar@1:       xfree(w);
alpar@1:       xfree(sol);
alpar@1:       xfree(vec);
alpar@1:       /* return to the calling program */
alpar@1:       return len;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: -- lpx_delete_cog - delete the conflict graph.
alpar@1: --
alpar@1: -- SYNOPSIS
alpar@1: --
alpar@1: -- #include "glplpx.h"
alpar@1: -- void lpx_delete_cog(void *cog);
alpar@1: --
alpar@1: -- DESCRIPTION
alpar@1: --
alpar@1: -- The routine lpx_delete_cog deletes the conflict graph, which the
alpar@1: -- parameter cog points to, freeing all the memory allocated to this
alpar@1: -- object. */
alpar@1: 
alpar@1: static void lpx_delete_cog(void *_cog)
alpar@1: {     struct COG *cog = _cog;
alpar@1:       xfree(cog->vert);
alpar@1:       xfree(cog->orig);
alpar@1:       xfree(cog->a);
alpar@1:       xfree(cog);
alpar@1: }
alpar@1: 
alpar@1: /**********************************************************************/
alpar@1: 
alpar@1: void *ios_clq_init(glp_tree *tree)
alpar@1: {     /* initialize clique cut generator */
alpar@1:       glp_prob *mip = tree->mip;
alpar@1:       xassert(mip != NULL);
alpar@1:       return lpx_create_cog(mip);
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  ios_clq_gen - generate clique cuts
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpios.h"
alpar@1: *  void ios_clq_gen(glp_tree *tree, void *gen);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine ios_clq_gen generates clique cuts for the current point
alpar@1: *  and adds them to the clique pool. */
alpar@1: 
alpar@1: void ios_clq_gen(glp_tree *tree, void *gen)
alpar@1: {     int n = lpx_get_num_cols(tree->mip);
alpar@1:       int len, *ind;
alpar@1:       double *val;
alpar@1:       xassert(gen != NULL);
alpar@1:       ind = xcalloc(1+n, sizeof(int));
alpar@1:       val = xcalloc(1+n, sizeof(double));
alpar@1:       len = lpx_clique_cut(tree->mip, gen, ind, val);
alpar@1:       if (len > 0)
alpar@1:       {  /* xprintf("len = %d\n", len); */
alpar@1:          glp_ios_add_row(tree, NULL, GLP_RF_CLQ, 0, len, ind, val,
alpar@1:             GLP_UP, val[0]);
alpar@1:       }
alpar@1:       xfree(ind);
alpar@1:       xfree(val);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /**********************************************************************/
alpar@1: 
alpar@1: void ios_clq_term(void *gen)
alpar@1: {     /* terminate clique cut generator */
alpar@1:       xassert(gen != NULL);
alpar@1:       lpx_delete_cog(gen);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /* eof */