/* cfg1.c (conflict graph) */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2012-2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * * 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 . ***********************************************************************/ #include "cfg.h" #include "env.h" #include "prob.h" #include "wclique.h" #include "wclique1.h" /*********************************************************************** * cfg_build_graph - build conflict graph * * This routine builds the conflict graph. It analyzes the specified * problem object to discover original and implied packing inequalities * and adds corresponding cliques to the conflict graph. * * Packing inequality has the form: * * sum z[j] <= 1, (1) * j in J * * where z[j] = x[j] or z[j] = 1 - x[j], x[j] is an original binary * variable. Every packing inequality (1) is equivalent to a set of * edge inequalities: * * z[i] + z[j] <= 1 for all i, j in J, i != j, (2) * * and since every edge inequality (2) defines an edge in the conflict * graph, corresponding packing inequality (1) defines a clique. * * To discover packing inequalities the routine analyzes constraints * of the specified MIP. To simplify the analysis each constraint is * analyzed separately. The analysis is performed as follows. * * Let some original constraint be the following: * * L <= sum a[j] x[j] <= U. (3) * * To analyze it the routine analyzes two constraints of "not greater * than" type: * * sum (-a[j]) x[j] <= -L, (4) * * sum (+a[j]) x[j] <= +U, (5) * * which are relaxations of the original constraint (3). (If, however, * L = -oo, or U = +oo, corresponding constraint being redundant is not * analyzed.) * * Let a constraint of "not greater than" type be the following: * * sum a[j] x[j] + sum a[j] x[j] <= b, (6) * j in J j in J' * * where J is a subset of binary variables, J' is a subset of other * (continues and non-binary integer) variables. The constraint (6) is * is relaxed as follows, to eliminate non-binary variables: * * sum a[j] x[j] <= b - sum a[j] x[j] <= b', (7) * j in J j in J' * * b' = sup(b - sum a[j] x[j]) = * j in J' * * = b - inf(sum a[j] x[j]) = * * = b - sum inf(a[j] x[j]) = (8) * * = b - sum a[j] inf(x[j]) - sum a[j] sup(x[j]) = * a[j]>0 a[j]<0 * * = b - sum a[j] l[j] - sum a[j] u[j], * a[j]>0 a[j]<0 * * where l[j] and u[j] are, resp., lower and upper bounds of x[j]. * * Then the routine transforms the relaxed constraint containing only * binary variables: * * sum a[j] x[j] <= b (9) * * to an equivalent 0-1 knapsack constraint as follows: * * sum a[j] x[j] + sum a[j] x[j] <= b ==> * a[j]>0 a[j]<0 * * sum a[j] x[j] + sum a[j] (1 - x[j]) <= b ==> * a[j]>0 a[j]<0 (10) * * sum (+a[j]) x[j] + sum (-a[j]) x[j] <= b + sum (-a[j]) ==> * a[j]>0 a[j]<0 a[j]<0 * * sum a'[j] z[j] <= b', * * where a'[j] = |a[j]| > 0, and * * ( x[j] if a[j] > 0 * z[j] = < * ( 1 - x[j] if a[j] < 0 * * is a binary variable, which is either original binary variable x[j] * or its complement. * * Finally, the routine analyzes the resultant 0-1 knapsack inequality: * * sum a[j] z[j] <= b, (11) * j in J * * where all a[j] are positive, to discover clique inequalities (1), * which are valid for (11) and therefore valid for (3). (It is assumed * that the original MIP has been preprocessed, so it is not checked, * for example, that b > 0 or that a[j] <= b.) * * In principle, to discover any edge inequalities valid for (11) it * is sufficient to check whether a[i] + a[j] > b for all i, j in J, * i < j. However, this way requires O(|J|^2) checks, so the routine * analyses (11) in the following way, which is much more efficient in * many practical cases. * * 1. Let a[p] and a[q] be two minimal coefficients: * * a[p] = min a[j], (12) * * a[q] = min a[j], j != p, (13) * * such that * * a[p] + a[q] > b. (14) * * This means that a[i] + a[j] > b for any i, j in J, i != j, so * * z[i] + z[j] <= 1 (15) * * are valid for (11) for any i, j in J, i != j. This case means that * J define a clique in the conflict graph. * * 2. Otherwise, let a[p] and [q] be two maximal coefficients: * * a[p] = max a[j], (16) * * a[q] = max a[j], j != p, (17) * * such that * * a[p] + a[q] <= b. (18) * * This means that a[i] + a[j] <= b for any i, j in J, i != j, so in * this case no valid edge inequalities for (11) exist. * * 3. Otherwise, let all a[j] be ordered by descending their values: * * a[1] >= a[2] >= ... >= a[p-1] >= a[p] >= a[p+1] >= ... (19) * * where p is such that * * a[p-1] + a[p] > b, (20) * * a[p] + a[p+1] <= b. (21) * * (May note that due to the former two cases in this case we always * have 2 <= p <= |J|-1.) * * Since a[p] and a[p-1] are two minimal coefficients in the set * J' = {1, ..., p}, J' define a clique in the conflict graph for the * same reason as in the first case. Similarly, since a[p] and a[p+1] * are two maximal coefficients in the set J" = {p, ..., |J|}, no edge * inequalities exist for all i, j in J" for the same reason as in the * second case. Thus, to discover other edge inequalities (15) valid * for (11), the routine checks if a[i] + a[j] > b for all i in J', * j in J", i != j. */ #define is_binary(j) \ (P->col[j]->kind == GLP_IV && P->col[j]->type == GLP_DB && \ P->col[j]->lb == 0.0 && P->col[j]->ub == 1.0) /* check if x[j] is binary variable */ struct term { int ind; double val; }; /* term a[j] * z[j] used to sort a[j]'s */ static int CDECL fcmp(const void *e1, const void *e2) { /* auxiliary routine called from qsort */ const struct term *t1 = e1, *t2 = e2; if (t1->val > t2->val) return -1; else if (t1->val < t2->val) return +1; else return 0; } static void analyze_ineq(glp_prob *P, CFG *G, int len, int ind[], double val[], double rhs, struct term t[]) { /* analyze inequality constraint (6) */ /* P is the original MIP * G is the conflict graph to be built * len is the number of terms in the constraint * ind[1], ..., ind[len] are indices of variables x[j] * val[1], ..., val[len] are constraint coefficients a[j] * rhs is the right-hand side b * t[1+len] is a working array */ int j, k, kk, p, q, type, new_len; /* eliminate non-binary variables; see (7) and (8) */ new_len = 0; for (k = 1; k <= len; k++) { /* get index of variable x[j] */ j = ind[k]; if (is_binary(j)) { /* x[j] remains in relaxed constraint */ new_len++; ind[new_len] = j; val[new_len] = val[k]; } else if (val[k] > 0.0) { /* eliminate non-binary x[j] in case a[j] > 0 */ /* b := b - a[j] * l[j]; see (8) */ type = P->col[j]->type; if (type == GLP_FR || type == GLP_UP) { /* x[j] has no lower bound */ goto done; } rhs -= val[k] * P->col[j]->lb; } else /* val[j] < 0.0 */ { /* eliminate non-binary x[j] in case a[j] < 0 */ /* b := b - a[j] * u[j]; see (8) */ type = P->col[j]->type; if (type == GLP_FR || type == GLP_LO) { /* x[j] has no upper bound */ goto done; } rhs -= val[k] * P->col[j]->ub; } } len = new_len; /* now we have the constraint (9) */ if (len <= 1) { /* at least two terms are needed */ goto done; } /* make all constraint coefficients positive; see (10) */ for (k = 1; k <= len; k++) { if (val[k] < 0.0) { /* a[j] < 0; substitute x[j] = 1 - x'[j], where x'[j] is * a complement binary variable */ ind[k] = -ind[k]; val[k] = -val[k]; rhs += val[k]; } } /* now we have 0-1 knapsack inequality (11) */ /* increase the right-hand side a bit to avoid false checks due * to rounding errors */ rhs += 0.001 * (1.0 + fabs(rhs)); /*** first case ***/ /* find two minimal coefficients a[p] and a[q] */ p = 0; for (k = 1; k <= len; k++) { if (p == 0 || val[p] > val[k]) p = k; } q = 0; for (k = 1; k <= len; k++) { if (k != p && (q == 0 || val[q] > val[k])) q = k; } xassert(p != 0 && q != 0 && p != q); /* check condition (14) */ if (val[p] + val[q] > rhs) { /* all z[j] define a clique in the conflict graph */ cfg_add_clique(G, len, ind); goto done; } /*** second case ***/ /* find two maximal coefficients a[p] and a[q] */ p = 0; for (k = 1; k <= len; k++) { if (p == 0 || val[p] < val[k]) p = k; } q = 0; for (k = 1; k <= len; k++) { if (k != p && (q == 0 || val[q] < val[k])) q = k; } xassert(p != 0 && q != 0 && p != q); /* check condition (18) */ if (val[p] + val[q] <= rhs) { /* no valid edge inequalities exist */ goto done; } /*** third case ***/ xassert(len >= 3); /* sort terms in descending order of coefficient values */ for (k = 1; k <= len; k++) { t[k].ind = ind[k]; t[k].val = val[k]; } qsort(&t[1], len, sizeof(struct term), fcmp); for (k = 1; k <= len; k++) { ind[k] = t[k].ind; val[k] = t[k].val; } /* now a[1] >= a[2] >= ... >= a[len-1] >= a[len] */ /* note that a[1] + a[2] > b and a[len-1] + a[len] <= b due two * the former two cases */ xassert(val[1] + val[2] > rhs); xassert(val[len-1] + val[len] <= rhs); /* find p according to conditions (20) and (21) */ for (p = 2; p < len; p++) { if (val[p] + val[p+1] <= rhs) break; } xassert(p < len); /* z[1], ..., z[p] define a clique in the conflict graph */ cfg_add_clique(G, p, ind); /* discover other edge inequalities */ for (k = 1; k <= p; k++) { for (kk = p; kk <= len; kk++) { if (k != kk && val[k] + val[kk] > rhs) { int iii[1+2]; iii[1] = ind[k]; iii[2] = ind[kk]; cfg_add_clique(G, 2, iii); } } } done: return; } CFG *cfg_build_graph(void *P_) { glp_prob *P = P_; int m = P->m; int n = P->n; CFG *G; int i, k, type, len, *ind; double *val; struct term *t; /* create the conflict graph (number of its vertices cannot be * greater than double number of binary variables) */ G = cfg_create_graph(n, 2 * glp_get_num_bin(P)); /* allocate working arrays */ ind = talloc(1+n, int); val = talloc(1+n, double); t = talloc(1+n, struct term); /* analyze constraints to discover edge inequalities */ for (i = 1; i <= m; i++) { type = P->row[i]->type; if (type == GLP_LO || type == GLP_DB || type == GLP_FX) { /* i-th row has lower bound */ /* analyze inequality sum (-a[j]) * x[j] <= -lb */ len = glp_get_mat_row(P, i, ind, val); for (k = 1; k <= len; k++) val[k] = -val[k]; analyze_ineq(P, G, len, ind, val, -P->row[i]->lb, t); } if (type == GLP_UP || type == GLP_DB || type == GLP_FX) { /* i-th row has upper bound */ /* analyze inequality sum (+a[j]) * x[j] <= +ub */ len = glp_get_mat_row(P, i, ind, val); analyze_ineq(P, G, len, ind, val, +P->row[i]->ub, t); } } /* free working arrays */ tfree(ind); tfree(val); tfree(t); return G; } /*********************************************************************** * cfg_find_clique - find maximum weight clique in conflict graph * * This routine finds a maximum weight clique in the conflict graph * G = (V, E), where the weight of vertex v in V is the value of * corresponding binary variable z (which is either an original binary * variable or its complement) in the optimal solution to LP relaxation * provided in the problem object. The goal is to find a clique in G, * whose weight is greater than 1, in which case corresponding packing * inequality is violated at the optimal point. * * On exit the routine stores vertex indices of the conflict graph * included in the clique found to locations ind[1], ..., ind[len], and * returns len, which is the clique size. The clique weight is stored * in location pointed to by the parameter sum. If no clique has been * found, the routine returns 0. * * Since the conflict graph may have a big number of vertices and be * quite dense, the routine uses an induced subgraph G' = (V', E'), * which is constructed as follows: * * 1. If the weight of some vertex v in V is zero (close to zero), it * is not included in V'. Obviously, including in a clique * zero-weight vertices does not change its weight, so if in G there * exist a clique of a non-zero weight, in G' exists a clique of the * same weight. This point is extremely important, because dropping * out zero-weight vertices can be done without retrieving lists of * adjacent vertices whose size may be very large. * * 2. Cumulative weight of vertex v in V is the sum of the weight of v * and weights of all vertices in V adjacent to v. Obviously, if * a clique includes a vertex v, the clique weight cannot be greater * than the cumulative weight of v. Since we are interested only in * cliques whose weight is greater than 1, vertices of V, whose * cumulative weight is not greater than 1, are not included in V'. * * May note that in many practical cases the size of the induced * subgraph G' is much less than the size of the original conflict * graph G due to many binary variables, whose optimal values are zero * or close to zero. For example, it may happen that |V| = 100,000 and * |E| = 1e9 while |V'| = 50 and |E'| = 1000. */ struct csa { /* common storage area */ glp_prob *P; /* original MIP */ CFG *G; /* original conflict graph G = (V, E), |V| = nv */ int *ind; /* int ind[1+nv]; */ /* working array */ /*--------------------------------------------------------------*/ /* induced subgraph G' = (V', E') of original conflict graph */ int nn; /* number of vertices in V' */ int *vtoi; /* int vtoi[1+nv]; */ /* vtoi[v] = i, 1 <= v <= nv, means that vertex v in V is vertex * i in V'; vtoi[v] = 0 means that vertex v is not included in * the subgraph */ int *itov; /* int itov[1+nv]; */ /* itov[i] = v, 1 <= i <= nn, means that vertex i in V' is vertex * v in V */ double *wgt; /* double wgt[1+nv]; */ /* wgt[i], 1 <= i <= nn, is a weight of vertex i in V', which is * the value of corresponding binary variable in optimal solution * to LP relaxation */ }; static void build_subgraph(struct csa *csa) { /* build induced subgraph */ glp_prob *P = csa->P; int n = P->n; CFG *G = csa->G; int *ind = csa->ind; int *pos = G->pos; int *neg = G->neg; int nv = G->nv; int *ref = G->ref; int *vtoi = csa->vtoi; int *itov = csa->itov; double *wgt = csa->wgt; int j, k, v, w, nn, len; double z, sum; /* initially induced subgraph is empty */ nn = 0; /* walk thru vertices of original conflict graph */ for (v = 1; v <= nv; v++) { /* determine value of binary variable z[j] that corresponds to * vertex v */ j = ref[v]; xassert(1 <= j && j <= n); if (pos[j] == v) { /* z[j] = x[j], where x[j] is original variable */ z = P->col[j]->prim; } else if (neg[j] == v) { /* z[j] = 1 - x[j], where x[j] is original variable */ z = 1.0 - P->col[j]->prim; } else xassert(v != v); /* if z[j] is close to zero, do not include v in the induced * subgraph */ if (z < 0.001) { vtoi[v] = 0; continue; } /* calculate cumulative weight of vertex v */ sum = z; /* walk thru all vertices adjacent to v */ len = cfg_get_adjacent(G, v, ind); for (k = 1; k <= len; k++) { /* there is an edge (v,w) in the conflict graph */ w = ind[k]; xassert(w != v); /* add value of z[j] that corresponds to vertex w */ j = ref[w]; xassert(1 <= j && j <= n); if (pos[j] == w) sum += P->col[j]->prim; else if (neg[j] == w) sum += 1.0 - P->col[j]->prim; else xassert(w != w); } /* cumulative weight of vertex v is an upper bound of weight * of any clique containing v; so if it not greater than 1, do * not include v in the induced subgraph */ if (sum < 1.010) { vtoi[v] = 0; continue; } /* include vertex v in the induced subgraph */ nn++; vtoi[v] = nn; itov[nn] = v; wgt[nn] = z; } /* induced subgraph has been built */ csa->nn = nn; return; } static int sub_adjacent(struct csa *csa, int i, int adj[]) { /* retrieve vertices of induced subgraph adjacent to specified * vertex */ CFG *G = csa->G; int nv = G->nv; int *ind = csa->ind; int nn = csa->nn; int *vtoi = csa->vtoi; int *itov = csa->itov; int j, k, v, w, len, len1; /* determine original vertex v corresponding to vertex i */ xassert(1 <= i && i <= nn); v = itov[i]; /* retrieve vertices adjacent to vertex v in original graph */ len1 = cfg_get_adjacent(G, v, ind); /* keep only adjacent vertices which are in induced subgraph and * change their numbers appropriately */ len = 0; for (k = 1; k <= len1; k++) { /* there exists edge (v, w) in original graph */ w = ind[k]; xassert(1 <= w && w <= nv && w != v); j = vtoi[w]; if (j != 0) { /* vertex w is vertex j in induced subgraph */ xassert(1 <= j && j <= nn && j != i); adj[++len] = j; } } return len; } static int find_clique(struct csa *csa, int c_ind[]) { /* find maximum weight clique in induced subgraph with exact * Ostergard's algorithm */ int nn = csa->nn; double *wgt = csa->wgt; int i, j, k, p, q, t, ne, nb, len, *iwt, *ind; unsigned char *a; xassert(nn >= 2); /* allocate working array */ ind = talloc(1+nn, int); /* calculate the number of elements in lower triangle (without * diagonal) of adjacency matrix of induced subgraph */ ne = (nn * (nn - 1)) / 2; /* calculate the number of bytes needed to store lower triangle * of adjacency matrix */ nb = (ne + (CHAR_BIT - 1)) / CHAR_BIT; /* allocate lower triangle of adjacency matrix */ a = talloc(nb, unsigned char); /* fill lower triangle of adjacency matrix */ memset(a, 0, nb); for (p = 1; p <= nn; p++) { /* retrieve vertices adjacent to vertex p */ len = sub_adjacent(csa, p, ind); for (k = 1; k <= len; k++) { /* there exists edge (p, q) in induced subgraph */ q = ind[k]; xassert(1 <= q && q <= nn && q != p); /* determine row and column indices of this edge in lower * triangle of adjacency matrix */ if (p > q) i = p, j = q; else /* p < q */ i = q, j = p; /* set bit a[i,j] to 1, i > j */ t = ((i - 1) * (i - 2)) / 2 + (j - 1); a[t / CHAR_BIT] |= (unsigned char)(1 << ((CHAR_BIT - 1) - t % CHAR_BIT)); } } /* scale vertex weights by 1000 and convert them to integers as * required by Ostergard's algorithm */ iwt = ind; for (i = 1; i <= nn; i++) { /* it is assumed that 0 <= wgt[i] <= 1 */ t = (int)(1000.0 * wgt[i] + 0.5); if (t < 0) t = 0; else if (t > 1000) t = 1000; iwt[i] = t; } /* find maximum weight clique */ len = wclique(nn, iwt, a, c_ind); /* free working arrays */ tfree(ind); tfree(a); /* return clique size to calling routine */ return len; } static int func(void *info, int i, int ind[]) { /* auxiliary routine used by routine find_clique1 */ struct csa *csa = info; xassert(1 <= i && i <= csa->nn); return sub_adjacent(csa, i, ind); } static int find_clique1(struct csa *csa, int c_ind[]) { /* find maximum weight clique in induced subgraph with greedy * heuristic */ int nn = csa->nn; double *wgt = csa->wgt; int len; xassert(nn >= 2); len = wclique1(nn, wgt, func, csa, c_ind); /* return clique size to calling routine */ return len; } int cfg_find_clique(void *P, CFG *G, int ind[], double *sum_) { int nv = G->nv; struct csa csa; int i, k, len; double sum; /* initialize common storage area */ csa.P = P; csa.G = G; csa.ind = talloc(1+nv, int); csa.nn = -1; csa.vtoi = talloc(1+nv, int); csa.itov = talloc(1+nv, int); csa.wgt = talloc(1+nv, double); /* build induced subgraph */ build_subgraph(&csa); #ifdef GLP_DEBUG xprintf("nn = %d\n", csa.nn); #endif /* if subgraph has less than two vertices, do nothing */ if (csa.nn < 2) { len = 0; sum = 0.0; goto skip; } /* find maximum weight clique in induced subgraph */ #if 1 /* FIXME */ if (csa.nn <= 50) #endif { /* induced subgraph is small; use exact algorithm */ len = find_clique(&csa, ind); } else { /* induced subgraph is large; use greedy heuristic */ len = find_clique1(&csa, ind); } /* do not report clique, if it has less than two vertices */ if (len < 2) { len = 0; sum = 0.0; goto skip; } /* convert indices of clique vertices from induced subgraph to * original conflict graph and compute clique weight */ sum = 0.0; for (k = 1; k <= len; k++) { i = ind[k]; xassert(1 <= i && i <= csa.nn); sum += csa.wgt[i]; ind[k] = csa.itov[i]; } skip: /* free working arrays */ tfree(csa.ind); tfree(csa.vtoi); tfree(csa.itov); tfree(csa.wgt); /* return to calling routine */ *sum_ = sum; return len; } /* eof */