diff options
Diffstat (limited to 'test/monniaux/glpk-4.65/src/intopt/cfg1.c')
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/cfg1.c | 703 |
1 files changed, 703 insertions, 0 deletions
diff --git a/test/monniaux/glpk-4.65/src/intopt/cfg1.c b/test/monniaux/glpk-4.65/src/intopt/cfg1.c new file mode 100644 index 00000000..80a2e834 --- /dev/null +++ b/test/monniaux/glpk-4.65/src/intopt/cfg1.c @@ -0,0 +1,703 @@ +/* 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: <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 "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 */ |