diff options
Diffstat (limited to 'test/monniaux/glpk-4.65/src/intopt')
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/cfg.c | 409 | ||||
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/cfg.h | 138 | ||||
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/cfg1.c | 703 | ||||
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/cfg2.c | 91 | ||||
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/clqcut.c | 134 | ||||
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/covgen.c | 885 | ||||
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/fpump.c | 360 | ||||
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/gmicut.c | 284 | ||||
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/gmigen.c | 142 | ||||
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/mirgen.c | 1529 | ||||
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/spv.c | 303 | ||||
| -rw-r--r-- | test/monniaux/glpk-4.65/src/intopt/spv.h | 83 |
12 files changed, 5061 insertions, 0 deletions
diff --git a/test/monniaux/glpk-4.65/src/intopt/cfg.c b/test/monniaux/glpk-4.65/src/intopt/cfg.c new file mode 100644 index 00000000..ab73b2da --- /dev/null +++ b/test/monniaux/glpk-4.65/src/intopt/cfg.c @@ -0,0 +1,409 @@ +/* cfg.c (conflict graph) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2012-2013 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" + +/*********************************************************************** +* cfg_create_graph - create conflict graph +* +* This routine creates the conflict graph, which initially is empty, +* and returns a pointer to the graph descriptor. +* +* The parameter n specifies the number of *all* variables in MIP, for +* which the conflict graph will be built. +* +* The parameter nv_max specifies maximal number of vertices in the +* conflict graph. It should be the double number of binary variables +* in corresponding MIP. */ + +CFG *cfg_create_graph(int n, int nv_max) +{ CFG *G; + xassert(n >= 0); + xassert(0 <= nv_max && nv_max <= n + n); + G = talloc(1, CFG); + G->n = n; + G->pos = talloc(1+n, int); + memset(&G->pos[1], 0, n * sizeof(int)); + G->neg = talloc(1+n, int); + memset(&G->neg[1], 0, n * sizeof(int)); + G->pool = dmp_create_pool(); + G->nv_max = nv_max; + G->nv = 0; + G->ref = talloc(1+nv_max, int); + G->vptr = talloc(1+nv_max, CFGVLE *); + G->cptr = talloc(1+nv_max, CFGCLE *); + return G; +} + +/*********************************************************************** +* cfg_add_clique - add clique to conflict graph +* +* This routine adds a clique to the conflict graph. +* +* The parameter size specifies the clique size, size >= 2. Note that +* any edge can be considered as a clique of size 2. +* +* The array ind specifies vertices constituting the clique in elements +* ind[k], 1 <= k <= size: +* +* ind[k] = +j means a vertex of the conflict graph that corresponds to +* original binary variable x[j], 1 <= j <= n. +* +* ind[k] = -j means a vertex of the conflict graph that corresponds to +* complement of original binary variable x[j], 1 <= j <= n. +* +* Note that if both vertices for x[j] and (1 - x[j]) have appeared in +* the conflict graph, the routine automatically adds an edge incident +* to these vertices. */ + +static void add_edge(CFG *G, int v, int w) +{ /* add clique of size 2 */ + DMP *pool = G->pool; + int nv = G->nv; + CFGVLE **vptr = G->vptr; + CFGVLE *vle; + xassert(1 <= v && v <= nv); + xassert(1 <= w && w <= nv); + xassert(v != w); + vle = dmp_talloc(pool, CFGVLE); + vle->v = w; + vle->next = vptr[v]; + vptr[v] = vle; + vle = dmp_talloc(pool, CFGVLE); + vle->v = v; + vle->next = vptr[w]; + vptr[w] = vle; + return; +} + +void cfg_add_clique(CFG *G, int size, const int ind[]) +{ int n = G->n; + int *pos = G->pos; + int *neg = G->neg; + DMP *pool = G->pool; + int nv_max = G->nv_max; + int *ref = G->ref; + CFGVLE **vptr = G->vptr; + CFGCLE **cptr = G->cptr; + int j, k, v; + xassert(2 <= size && size <= nv_max); + /* add new vertices to the conflict graph */ + for (k = 1; k <= size; k++) + { j = ind[k]; + if (j > 0) + { /* vertex corresponds to x[j] */ + xassert(1 <= j && j <= n); + if (pos[j] == 0) + { /* no such vertex exists; add it */ + v = pos[j] = ++(G->nv); + xassert(v <= nv_max); + ref[v] = j; + vptr[v] = NULL; + cptr[v] = NULL; + if (neg[j] != 0) + { /* now both vertices for x[j] and (1 - x[j]) exist */ + add_edge(G, v, neg[j]); + } + } + } + else + { /* vertex corresponds to (1 - x[j]) */ + j = -j; + xassert(1 <= j && j <= n); + if (neg[j] == 0) + { /* no such vertex exists; add it */ + v = neg[j] = ++(G->nv); + xassert(v <= nv_max); + ref[v] = j; + vptr[v] = NULL; + cptr[v] = NULL; + if (pos[j] != 0) + { /* now both vertices for x[j] and (1 - x[j]) exist */ + add_edge(G, v, pos[j]); + } + } + } + } + /* add specified clique to the conflict graph */ + if (size == 2) + add_edge(G, + ind[1] > 0 ? pos[+ind[1]] : neg[-ind[1]], + ind[2] > 0 ? pos[+ind[2]] : neg[-ind[2]]); + else + { CFGVLE *vp, *vle; + CFGCLE *cle; + /* build list of clique vertices */ + vp = NULL; + for (k = 1; k <= size; k++) + { vle = dmp_talloc(pool, CFGVLE); + vle->v = ind[k] > 0 ? pos[+ind[k]] : neg[-ind[k]]; + vle->next = vp; + vp = vle; + } + /* attach the clique to all its vertices */ + for (k = 1; k <= size; k++) + { cle = dmp_talloc(pool, CFGCLE); + cle->vptr = vp; + v = ind[k] > 0 ? pos[+ind[k]] : neg[-ind[k]]; + cle->next = cptr[v]; + cptr[v] = cle; + } + } + return; +} + +/*********************************************************************** +* cfg_get_adjacent - get vertices adjacent to specified vertex +* +* This routine stores numbers of all vertices adjacent to specified +* vertex v of the conflict graph in locations ind[1], ..., ind[len], +* and returns len, 1 <= len <= nv-1, where nv is the total number of +* vertices in the conflict graph. +* +* Note that the conflict graph defined by this routine has neither +* self-loops nor multiple edges. */ + +int cfg_get_adjacent(CFG *G, int v, int ind[]) +{ int nv = G->nv; + int *ref = G->ref; + CFGVLE **vptr = G->vptr; + CFGCLE **cptr = G->cptr; + CFGVLE *vle; + CFGCLE *cle; + int k, w, len; + xassert(1 <= v && v <= nv); + len = 0; + /* walk thru the list of adjacent vertices */ + for (vle = vptr[v]; vle != NULL; vle = vle->next) + { w = vle->v; + xassert(1 <= w && w <= nv); + xassert(w != v); + if (ref[w] > 0) + { ind[++len] = w; + ref[w] = -ref[w]; + } + } + /* walk thru the list of incident cliques */ + for (cle = cptr[v]; cle != NULL; cle = cle->next) + { /* walk thru the list of clique vertices */ + for (vle = cle->vptr; vle != NULL; vle = vle->next) + { w = vle->v; + xassert(1 <= w && w <= nv); + if (w != v && ref[w] > 0) + { ind[++len] = w; + ref[w] = -ref[w]; + } + } + } + xassert(1 <= len && len < nv); + /* unmark vertices included in the resultant adjacency list */ + for (k = 1; k <= len; k++) + { w = ind[k]; + ref[w] = -ref[w]; + } + return len; +} + +/*********************************************************************** +* cfg_expand_clique - expand specified clique to maximal clique +* +* Given some clique in the conflict graph this routine expands it to +* a maximal clique by including in it new vertices. +* +* On entry vertex indices constituting the initial clique should be +* stored in locations c_ind[1], ..., c_ind[c_len], where c_len is the +* initial clique size. On exit the routine stores new vertex indices +* to locations c_ind[c_len+1], ..., c_ind[c_len'], where c_len' is the +* size of the maximal clique found, and returns c_len'. +* +* ALGORITHM +* +* Let G = (V, E) be a graph, C within V be a current clique to be +* expanded, and D within V \ C be a subset of vertices adjacent to all +* vertices from C. On every iteration the routine chooses some vertex +* v in D, includes it into C, and removes from D the vertex v as well +* as all vertices not adjacent to v. Initially C is empty and D = V. +* Iterations repeat until D becomes an empty set. Obviously, the final +* set C is a maximal clique in G. +* +* Now let C0 be an initial clique, and we want C0 to be a subset of +* the final maximal clique C. To provide this condition the routine +* starts constructing C by choosing only such vertices v in D, which +* are in C0, until all vertices from C0 have been included in C. May +* note that if on some iteration C0 \ C is non-empty (i.e. if not all +* vertices from C0 have been included in C), C0 \ C is a subset of D, +* because C0 is a clique. */ + +static int intersection(int d_len, int d_ind[], int d_pos[], int len, + const int ind[]) +{ /* compute intersection D := D inter W, where W is some specified + * set of vertices */ + int k, t, v, new_len; + /* walk thru vertices in W and mark vertices in D */ + for (t = 1; t <= len; t++) + { /* v in W */ + v = ind[t]; + /* determine position of v in D */ + k = d_pos[v]; + if (k != 0) + { /* v in D */ + xassert(d_ind[k] == v); + /* mark v to keep it in D */ + d_ind[k] = -v; + } + } + /* remove all unmarked vertices from D */ + new_len = 0; + for (k = 1; k <= d_len; k++) + { /* v in D */ + v = d_ind[k]; + if (v < 0) + { /* v is marked; keep it */ + v = -v; + new_len++; + d_ind[new_len] = v; + d_pos[v] = new_len; + } + else + { /* v is not marked; remove it */ + d_pos[v] = 0; + } + } + return new_len; +} + +int cfg_expand_clique(CFG *G, int c_len, int c_ind[]) +{ int nv = G->nv; + int d_len, *d_ind, *d_pos, len, *ind; + int k, v; + xassert(0 <= c_len && c_len <= nv); + /* allocate working arrays */ + d_ind = talloc(1+nv, int); + d_pos = talloc(1+nv, int); + ind = talloc(1+nv, int); + /* initialize C := 0, D := V */ + d_len = nv; + for (k = 1; k <= nv; k++) + d_ind[k] = d_pos[k] = k; + /* expand C by vertices of specified initial clique C0 */ + for (k = 1; k <= c_len; k++) + { /* v in C0 */ + v = c_ind[k]; + xassert(1 <= v && v <= nv); + /* since C0 is clique, v should be in D */ + xassert(d_pos[v] != 0); + /* W := set of vertices adjacent to v */ + len = cfg_get_adjacent(G, v, ind); + /* D := D inter W */ + d_len = intersection(d_len, d_ind, d_pos, len, ind); + /* since v not in W, now v should be not in D */ + xassert(d_pos[v] == 0); + } + /* expand C by some other vertices until D is empty */ + while (d_len > 0) + { /* v in D */ + v = d_ind[1]; + xassert(1 <= v && v <= nv); + /* note that v is adjacent to all vertices in C (by design), + * so add v to C */ + c_ind[++c_len] = v; + /* W := set of vertices adjacent to v */ + len = cfg_get_adjacent(G, v, ind); + /* D := D inter W */ + d_len = intersection(d_len, d_ind, d_pos, len, ind); + /* since v not in W, now v should be not in D */ + xassert(d_pos[v] == 0); + } + /* free working arrays */ + tfree(d_ind); + tfree(d_pos); + tfree(ind); + /* bring maximal clique to calling routine */ + return c_len; +} + +/*********************************************************************** +* cfg_check_clique - check clique in conflict graph +* +* This routine checks that vertices of the conflict graph specified +* in locations c_ind[1], ..., c_ind[c_len] constitute a clique. +* +* NOTE: for testing/debugging only. */ + +void cfg_check_clique(CFG *G, int c_len, const int c_ind[]) +{ int nv = G->nv; + int k, kk, v, w, len, *ind; + char *flag; + ind = talloc(1+nv, int); + flag = talloc(1+nv, char); + memset(&flag[1], 0, nv); + /* walk thru clique vertices */ + xassert(c_len >= 0); + for (k = 1; k <= c_len; k++) + { /* get clique vertex v */ + v = c_ind[k]; + xassert(1 <= v && v <= nv); + /* get vertices adjacent to vertex v */ + len = cfg_get_adjacent(G, v, ind); + for (kk = 1; kk <= len; kk++) + { w = ind[kk]; + xassert(1 <= w && w <= nv); + xassert(w != v); + flag[w] = 1; + } + /* check that all clique vertices other than v are adjacent + to v */ + for (kk = 1; kk <= c_len; kk++) + { w = c_ind[kk]; + xassert(1 <= w && w <= nv); + if (w != v) + xassert(flag[w]); + } + /* reset vertex flags */ + for (kk = 1; kk <= len; kk++) + flag[ind[kk]] = 0; + } + tfree(ind); + tfree(flag); + return; +} + +/*********************************************************************** +* cfg_delete_graph - delete conflict graph +* +* This routine deletes the conflict graph by freeing all the memory +* allocated to this program object. */ + +void cfg_delete_graph(CFG *G) +{ tfree(G->pos); + tfree(G->neg); + dmp_delete_pool(G->pool); + tfree(G->ref); + tfree(G->vptr); + tfree(G->cptr); + tfree(G); + return; +} + +/* eof */ diff --git a/test/monniaux/glpk-4.65/src/intopt/cfg.h b/test/monniaux/glpk-4.65/src/intopt/cfg.h new file mode 100644 index 00000000..d478f6c0 --- /dev/null +++ b/test/monniaux/glpk-4.65/src/intopt/cfg.h @@ -0,0 +1,138 @@ +/* cfg.h (conflict graph) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2012-2013 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/>. +***********************************************************************/ + +#ifndef CFG_H +#define CFG_H + +#include "dmp.h" + +/*********************************************************************** +* The structure CFG describes the conflict graph. +* +* Conflict graph is an undirected graph G = (V, E), where V is a set +* of vertices, E <= V x V is a set of edges. Each vertex v in V of the +* conflict graph corresponds to a binary variable z[v], which is +* either an original binary variable x[j] or its complement 1 - x[j]. +* Edge (v,w) in E means that z[v] and z[w] cannot take the value 1 at +* the same time, i.e. it defines an inequality z[v] + z[w] <= 1, which +* is assumed to be valid for original MIP. +* +* Since the conflict graph may be dense, it is stored as an union of +* its cliques rather than explicitly. */ + +#if 0 /* 08/III-2016 */ +typedef struct CFG CFG; +#else +typedef struct glp_cfg CFG; +#endif +typedef struct CFGVLE CFGVLE; +typedef struct CFGCLE CFGCLE; + +#if 0 /* 08/III-2016 */ +struct CFG +#else +struct glp_cfg +#endif +{ /* conflict graph descriptor */ + int n; + /* number of *all* variables (columns) in corresponding MIP */ + int *pos; /* int pos[1+n]; */ + /* pos[0] is not used; + * pos[j] = v, 1 <= j <= n, means that vertex v corresponds to + * original binary variable x[j], and pos[j] = 0 means that the + * conflict graph has no such vertex */ + int *neg; /* int neg[1+n]; */ + /* neg[0] is not used; + * neg[j] = v, 1 <= j <= n, means that vertex v corresponds to + * complement of original binary variable x[j], and neg[j] = 0 + * means that the conflict graph has no such vertex */ + DMP *pool; + /* memory pool to allocate elements of the conflict graph */ + int nv_max; + /* maximal number of vertices in the conflict graph */ + int nv; + /* current number of vertices in the conflict graph */ + int *ref; /* int ref[1+nv_max]; */ + /* ref[v] = j, 1 <= v <= nv, means that vertex v corresponds + * either to original binary variable x[j] or to its complement, + * i.e. either pos[j] = v or neg[j] = v */ + CFGVLE **vptr; /* CFGVLE *vptr[1+nv_max]; */ + /* vptr[v], 1 <= v <= nv, is an initial pointer to the list of + * vertices adjacent to vertex v */ + CFGCLE **cptr; /* CFGCLE *cptr[1+nv_max]; */ + /* cptr[v], 1 <= v <= nv, is an initial pointer to the list of + * cliques that contain vertex v */ +}; + +struct CFGVLE +{ /* vertex list element */ + int v; + /* vertex number, 1 <= v <= nv */ + CFGVLE *next; + /* pointer to next vertex list element */ +}; + +struct CFGCLE +{ /* clique list element */ + CFGVLE *vptr; + /* initial pointer to the list of clique vertices */ + CFGCLE *next; + /* pointer to next clique list element */ +}; + +#define cfg_create_graph _glp_cfg_create_graph +CFG *cfg_create_graph(int n, int nv_max); +/* create conflict graph */ + +#define cfg_add_clique _glp_cfg_add_clique +void cfg_add_clique(CFG *G, int size, const int ind[]); +/* add clique to conflict graph */ + +#define cfg_get_adjacent _glp_cfg_get_adjacent +int cfg_get_adjacent(CFG *G, int v, int ind[]); +/* get vertices adjacent to specified vertex */ + +#define cfg_expand_clique _glp_cfg_expand_clique +int cfg_expand_clique(CFG *G, int c_len, int c_ind[]); +/* expand specified clique to maximal clique */ + +#define cfg_check_clique _glp_cfg_check_clique +void cfg_check_clique(CFG *G, int c_len, const int c_ind[]); +/* check clique in conflict graph */ + +#define cfg_delete_graph _glp_cfg_delete_graph +void cfg_delete_graph(CFG *G); +/* delete conflict graph */ + +#define cfg_build_graph _glp_cfg_build_graph +CFG *cfg_build_graph(void /* glp_prob */ *P); +/* build conflict graph */ + +#define cfg_find_clique _glp_cfg_find_clique +int cfg_find_clique(void /* glp_prob */ *P, CFG *G, int ind[], + double *sum); +/* find maximum weight clique in conflict graph */ + +#endif + +/* eof */ 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 */ diff --git a/test/monniaux/glpk-4.65/src/intopt/cfg2.c b/test/monniaux/glpk-4.65/src/intopt/cfg2.c new file mode 100644 index 00000000..85c0705e --- /dev/null +++ b/test/monniaux/glpk-4.65/src/intopt/cfg2.c @@ -0,0 +1,91 @@ +/* cfg2.c (conflict graph) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2015-2016 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" + +/*********************************************************************** +* NAME +* +* glp_cfg_init - create and initialize conflict graph +* +* SYNOPSIS +* +* glp_cfg *glp_cfg_init(glp_prob *P); +* +* DESCRIPTION +* +* This routine creates and initializes the conflict graph for the +* specified problem object. +* +* RETURNS +* +* The routine returns a pointer to the conflict graph descriptor. +* However, if the conflict graph is empty (no conflicts have been +* found), the routine returns NULL. */ + +glp_cfg *glp_cfg_init(glp_prob *P) +{ glp_cfg *G; + int j, n1, n2; + xprintf("Constructing conflict graph...\n"); + G = cfg_build_graph(P); + n1 = n2 = 0; + for (j = 1; j <= P->n; j++) + { if (G->pos[j]) + n1 ++; + if (G->neg[j]) + n2++; + } + if (n1 == 0 && n2 == 0) + { xprintf("No conflicts found\n"); + cfg_delete_graph(G); + G = NULL; + } + else + xprintf("Conflict graph has %d + %d = %d vertices\n", + n1, n2, G->nv); + return G; +} + +/*********************************************************************** +* NAME +* +* glp_cfg_free - delete conflict graph descriptor +* +* SYNOPSIS +* +* void glp_cfg_free(glp_cfg *G); +* +* DESCRIPTION +* +* This routine deletes the conflict graph descriptor and frees all the +* memory allocated to it. */ + +void glp_cfg_free(glp_cfg *G) +{ xassert(G != NULL); + cfg_delete_graph(G); + return; +} + +/* eof */ diff --git a/test/monniaux/glpk-4.65/src/intopt/clqcut.c b/test/monniaux/glpk-4.65/src/intopt/clqcut.c new file mode 100644 index 00000000..d3db5b39 --- /dev/null +++ b/test/monniaux/glpk-4.65/src/intopt/clqcut.c @@ -0,0 +1,134 @@ +/* clqcut.c (clique cut generator) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2008-2016 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" + +/*********************************************************************** +* NAME +* +* glp_clq_cut - generate clique cut from conflict graph +* +* SYNOPSIS +* +* int glp_clq_cut(glp_prob *P, glp_cfg *G, int ind[], double val[]); +* +* DESCRIPTION +* +* This routine attempts to generate a clique cut. +* +* The cut generated by the routine is the following inequality: +* +* sum a[j] * x[j] <= b, +* +* which is expected to be violated at the current basic solution. +* +* If the cut has been successfully generated, the routine stores its +* non-zero coefficients a[j] and corresponding column indices j in the +* array locations val[1], ..., val[len] and ind[1], ..., ind[len], +* where 1 <= len <= n is the number of non-zero coefficients. The +* right-hand side value b is stored in val[0], and ind[0] is set to 0. +* +* RETURNS +* +* If the cut has been successfully generated, the routine returns +* len, the number of non-zero coefficients in the cut, 1 <= len <= n. +* Otherwise, the routine returns a non-positive value. */ + +int glp_clq_cut(glp_prob *P, glp_cfg *G, int ind[], double val[]) +{ int n = P->n; + int *pos = G->pos; + int *neg = G->neg; + int nv = G->nv; + int *ref = G->ref; + int j, k, v, len; + double rhs, sum; + xassert(G->n == n); + /* find maximum weight clique in conflict graph */ + len = cfg_find_clique(P, G, ind, &sum); +#ifdef GLP_DEBUG + xprintf("len = %d; sum = %g\n", len, sum); + cfg_check_clique(G, len, ind); +#endif + /* check if clique inequality is violated */ + if (sum < 1.07) + return 0; + /* expand clique to maximal one */ + len = cfg_expand_clique(G, len, ind); +#ifdef GLP_DEBUG + xprintf("maximal clique size = %d\n", len); + cfg_check_clique(G, len, ind); +#endif + /* construct clique cut (fixed binary variables are removed, so + this cut is only locally valid) */ + rhs = 1.0; + for (j = 1; j <= n; j++) + val[j] = 0.0; + for (k = 1; k <= len; k++) + { /* v is clique vertex */ + v = ind[k]; + xassert(1 <= v && v <= nv); + /* j is number of corresponding binary variable */ + j = ref[v]; + xassert(1 <= j && j <= n); + if (pos[j] == v) + { /* v corresponds to x[j] */ + if (P->col[j]->type == GLP_FX) + { /* x[j] is fixed */ + rhs -= P->col[j]->prim; + } + else + { /* x[j] is not fixed */ + val[j] += 1.0; + } + } + else if (neg[j] == v) + { /* v corresponds to (1 - x[j]) */ + if (P->col[j]->type == GLP_FX) + { /* x[j] is fixed */ + rhs -= (1.0 - P->col[j]->prim); + } + else + { /* x[j] is not fixed */ + val[j] -= 1.0; + rhs -= 1.0; + } + } + else + xassert(v != v); + } + /* convert cut inequality to sparse format */ + len = 0; + for (j = 1; j <= n; j++) + { if (val[j] != 0.0) + { len++; + ind[len] = j; + val[len] = val[j]; + } + } + ind[0] = 0, val[0] = rhs; + return len; +} + +/* eof */ diff --git a/test/monniaux/glpk-4.65/src/intopt/covgen.c b/test/monniaux/glpk-4.65/src/intopt/covgen.c new file mode 100644 index 00000000..427c3aa8 --- /dev/null +++ b/test/monniaux/glpk-4.65/src/intopt/covgen.c @@ -0,0 +1,885 @@ +/* covgen.c */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2017-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 "env.h" +#include "fvs.h" +#include "ks.h" +#include "prob.h" + +struct glp_cov +{ /* cover cut generator working area */ + int n; + /* number of columns (variables) */ + glp_prob *set; + /* set of globally valid 0-1 knapsack inequalities chosen from + * the root problem; each inequality is either original row or + * its relaxation (surrogate 0-1 knapsack) which is constructed + * by substitution of lower/upper single/variable bounds for + * continuous and general integer (non-binary) variables */ +}; + +struct bnd +{ /* simple or variable bound */ + /* if z = 0, it is a simple bound x >= or <= b; if b = -DBL_MAX + * (b = +DBL_MAX), x has no lower (upper) bound; otherwise, if + * z != 0, it is a variable bound x >= or <= a * z + b */ + int z; + /* number of binary variable or 0 */ + double a, b; + /* bound parameters */ +}; + +struct csa +{ /* common storage area */ + glp_prob *P; + /* original (root) MIP */ + struct bnd *l; /* struct bnd l[1+P->n]; */ + /* lower simple/variable bounds of variables */ + struct bnd *u; /* struct bnd u[1+P->n]; */ + /* upper simple/variable bounds of variables */ + glp_prob *set; + /* see struct glp_cov above */ +}; + +/*********************************************************************** +* init_bounds - initialize bounds of variables with simple bounds +* +* This routine initializes lower and upper bounds of all variables +* with simple bounds specified in the original mip. */ + +static void init_bounds(struct csa *csa) +{ glp_prob *P = csa->P; + struct bnd *l = csa->l, *u = csa->u; + int j; + for (j = 1; j <= P->n; j++) + { l[j].z = u[j].z = 0; + l[j].a = u[j].a = 0; + l[j].b = glp_get_col_lb(P, j); + u[j].b = glp_get_col_ub(P, j); + } + return; +} + +/*********************************************************************** +* check_vb - check variable bound +* +* This routine checks if the specified i-th row has the form +* +* a1 * x + a2 * z >= or <= rhs, (1) +* +* where x is a non-fixed continuous or general integer variable, and +* z is a binary variable. If it is, the routine converts the row to +* the following variable lower/upper bound (VLB/VUB) of x: +* +* x >= or <= a * z + b, (2) +* +* where a = - a2 / a1, b = rhs / a1. Note that the inequality type is +* changed to opposite one when a1 < 0. +* +* If the row is identified as a variable bound, the routine returns +* GLP_LO for VLB or GLP_UP for VUB and provides the reference numbers +* of variables x and z and values of a and b. Otherwise, the routine +* returns zero. */ + +static int check_vb(struct csa *csa, int i, int *x, int *z, double *a, + double *b) +{ glp_prob *P = csa->P; + GLPROW *row; + GLPAIJ *a1, *a2; + int type; + double rhs; + xassert(1 <= i && i <= P->m); + row = P->row[i]; + /* check row type */ + switch (row->type) + { case GLP_LO: + case GLP_UP: + break; + default: + return 0; + } + /* take first term of the row */ + a1 = row->ptr; + if (a1 == NULL) + return 0; + /* take second term of the row */ + a2 = a1->r_next; + if (a2 == NULL) + return 0; + /* there should be exactly two terms in the row */ + if (a2->r_next != NULL) + return 0; + /* if first term is a binary variable, swap the terms */ + if (glp_get_col_kind(P, a1->col->j) == GLP_BV) + { GLPAIJ *a; + a = a1, a1 = a2, a2 = a; + } + /* now first term should be a non-fixed continuous or general + * integer variable */ + if (a1->col->type == GLP_FX) + return 0; + if (glp_get_col_kind(P, a1->col->j) == GLP_BV) + return 0; + /* and second term should be a binary variable */ + if (glp_get_col_kind(P, a2->col->j) != GLP_BV) + return 0; + /* VLB/VUB row has been identified */ + switch (row->type) + { case GLP_LO: + type = a1->val > 0 ? GLP_LO : GLP_UP; + rhs = row->lb; + break; + case GLP_UP: + type = a1->val > 0 ? GLP_UP : GLP_LO; + rhs = row->ub; + break; + default: + xassert(type != type); + } + *x = a1->col->j; + *z = a2->col->j; + *a = - a2->val / a1->val; + *b = rhs / a1->val; + return type; +} + +/*********************************************************************** +* set_vb - set variable bound +* +* This routine sets lower or upper variable bound specified as +* +* x >= a * z + b (type = GLP_LO) +* +* x <= a * z + b (type = GLP_UP) */ + +static void set_vb(struct csa *csa, int type, int x, int z, double a, + double b) +{ glp_prob *P = csa->P; + struct bnd *l = csa->l, *u = csa->u; + xassert(glp_get_col_type(P, x) != GLP_FX); + xassert(glp_get_col_kind(P, x) != GLP_BV); + xassert(glp_get_col_kind(P, z) == GLP_BV); + xassert(a != 0); + switch (type) + { case GLP_LO: + /* FIXME: check existing simple lower bound? */ + l[x].z = z, l[x].a = a, l[x].b = b; + break; + case GLP_UP: + /* FIXME: check existing simple upper bound? */ + u[x].z = z, u[x].a = a, u[x].b = b; + break; + default: + xassert(type != type); + } + return; +} + +/*********************************************************************** +* obtain_vbs - obtain and set variable bounds +* +* This routine walks thru all rows of the original mip, identifies +* rows specifying variable lower/upper bounds, and sets these bounds +* for corresponding (non-binary) variables. */ + +static void obtain_vbs(struct csa *csa) +{ glp_prob *P = csa->P; + int i, x, z, type, save; + double a, b; + for (i = 1; i <= P->m; i++) + { switch (P->row[i]->type) + { case GLP_FR: + break; + case GLP_LO: + case GLP_UP: + type = check_vb(csa, i, &x, &z, &a, &b); + if (type) + set_vb(csa, type, x, z, a, b); + break; + case GLP_DB: + case GLP_FX: + /* double-side inequality l <= ... <= u and equality + * ... = l = u are considered as two single inequalities + * ... >= l and ... <= u */ + save = P->row[i]->type; + P->row[i]->type = GLP_LO; + type = check_vb(csa, i, &x, &z, &a, &b); + if (type) + set_vb(csa, type, x, z, a, b); + P->row[i]->type = GLP_UP; + type = check_vb(csa, i, &x, &z, &a, &b); + if (type) + set_vb(csa, type, x, z, a, b); + P->row[i]->type = save; + break; + default: + xassert(P != P); + } + } + return; +} + +/*********************************************************************** +* add_term - add term to sparse vector +* +* This routine computes the following linear combination: +* +* v := v + a * e[j], +* +* where v is a sparse vector in full storage format, a is a non-zero +* scalar, e[j] is j-th column of unity matrix. */ + +static void add_term(FVS *v, int j, double a) +{ xassert(1 <= j && j <= v->n); + xassert(a != 0); + if (v->vec[j] == 0) + { /* create j-th component */ + v->nnz++; + xassert(v->nnz <= v->n); + v->ind[v->nnz] = j; + } + /* perform addition */ + v->vec[j] += a; + if (fabs(v->vec[j]) < 1e-9 * (1 + fabs(a))) + { /* remove j-th component */ + v->vec[j] = DBL_MIN; + } + return; +} + +/*********************************************************************** +* build_ks - build "0-1 knapsack" inequality +* +* Given an inequality of "not greater" type: +* +* sum{j in 1..n} a[j]*x[j] <= b, (1) +* +* this routine attempts to transform it to equivalent or relaxed "0-1 +* knapsack" inequality that contains only binary variables. +* +* If x[j] is a binary variable, the term a[j]*x[j] is not changed. +* Otherwise, if x[j] is a continuous or integer non-binary variable, +* it is replaced by its lower (if a[j] > 0) or upper (if a[j] < 0) +* single or variable bound. In the latter case, if x[j] is a non-fixed +* variable, this results in a relaxation of original inequality known +* as "surrogate knapsack". Thus, if the specified inequality is valid +* for the original mip, the resulting inequality is also valid. +* +* Note that in both source and resulting inequalities coefficients +* a[j] can have any sign. +* +* On entry to the routine the source inequality is specified by the +* parameters n, ind (contains original numbers of x[j]), a, and b. The +* parameter v is a working sparse vector whose components are assumed +* to be zero. +* +* On exit the routine stores the resulting "0-1 knapsack" inequality +* in the parameters ind, a, and b, and returns n which is the number +* of terms in the resulting inequality. Zero content of the vector v +* is restored before exit. +* +* If the resulting inequality cannot be constructed due to missing +* lower/upper bounds of some variable, the routine returns a negative +* value. */ + +static int build_ks(struct csa *csa, int n, int ind[], double a[], + double *b, FVS *v) +{ glp_prob *P = csa->P; + struct bnd *l = csa->l, *u = csa->u; + int j, k; + /* check that v = 0 */ +#ifdef GLP_DEBUG + fvs_check_vec(v); +#endif + xassert(v->nnz == 0); + /* walk thru terms of original inequality */ + for (j = 1; j <= n; j++) + { /* process term a[j]*x[j] */ + k = ind[j]; /* original number of x[j] in mip */ + if (glp_get_col_kind(P, k) == GLP_BV) + { /* x[j] is a binary variable */ + /* include its term into resulting inequality */ + add_term(v, k, a[j]); + } + else if (a[j] > 0) + { /* substitute x[j] by its lower bound */ + if (l[k].b == -DBL_MAX) + { /* x[j] has no lower bound */ + n = -1; + goto skip; + } + else if (l[k].z == 0) + { /* x[j] has simple lower bound */ + *b -= a[j] * l[k].b; + } + else + { /* x[j] has variable lower bound (a * z + b) */ + add_term(v, l[k].z, a[j] * l[k].a); + *b -= a[j] * l[k].b; + } + } + else /* a[j] < 0 */ + { /* substitute x[j] by its upper bound */ + if (u[k].b == +DBL_MAX) + { /* x[j] has no upper bound */ + n = -1; + goto skip; + } + else if (u[k].z == 0) + { /* x[j] has simple upper bound */ + *b -= a[j] * u[k].b; + } + else + { /* x[j] has variable upper bound (a * z + b) */ + add_term(v, u[k].z, a[j] * u[k].a); + *b -= a[j] * u[k].b; + } + } + } + /* replace tiny coefficients by exact zeros (see add_term) */ + fvs_adjust_vec(v, 2 * DBL_MIN); + /* copy terms of resulting inequality */ + xassert(v->nnz <= n); + n = v->nnz; + for (j = 1; j <= n; j++) + { ind[j] = v->ind[j]; + a[j] = v->vec[ind[j]]; + } +skip: /* restore zero content of v */ + fvs_clear_vec(v); + return n; +} + +/*********************************************************************** +* can_be_active - check if inequality can be active +* +* This routine checks if the specified "0-1 knapsack" inequality +* +* sum{j in 1..n} a[j]*x[j] <= b +* +* can be active. If so, the routine returns true, otherwise false. */ + +static int can_be_active(int n, const double a[], double b) +{ int j; + double s; + s = 0; + for (j = 1; j <= n; j++) + { if (a[j] > 0) + s += a[j]; + } + return s > b + .001 * (1 + fabs(b)); +} + +/*********************************************************************** +* is_sos_ineq - check if inequality is packing (SOS) constraint +* +* This routine checks if the specified "0-1 knapsack" inequality +* +* sum{j in 1..n} a[j]*x[j] <= b (1) +* +* is equivalent to packing inequality (Padberg calls such inequalities +* special ordered set or SOS constraints) +* +* sum{j in J'} x[j] - sum{j in J"} x[j] <= 1 - |J"|. (2) +* +* If so, the routine returns true, otherwise false. +* +* Note that if X is a set of feasible binary points satisfying to (2), +* its convex hull conv(X) equals to the set of feasible points of LP +* relaxation of (2), which is a n-dimensional simplex, so inequalities +* (2) are useless for generating cover cuts (due to unimodularity). +* +* ALGORITHM +* +* First, we make all a[j] positive by complementing x[j] = 1 - x'[j] +* in (1). This is performed implicitly (i.e. actually the array a is +* not changed), but b is replaced by b - sum{j : a[j] < 0}. +* +* Then we find two smallest coefficients a[p] = min{j in 1..n} a[j] +* and a[q] = min{j in 1..n : j != p} a[j]. It is obvious that if +* a[p] + a[q] > b, then a[i] + a[j] > b for all i != j, from which it +* follows that x[i] + x[j] <= 1 for all i != j. But the latter means +* that the original inequality (with all a[j] > 0) is equivalent to +* packing inequality +* +* sum{j in 1..n} x[j] <= 1. (3) +* +* Returning to original (uncomplemented) variables x'[j] = 1 - x[j] +* we have that the original inequality is equivalent to (2), where +* J' = {j : a[j] > 0} and J" = {j : a[j] < 0}. */ + +static int is_sos_ineq(int n, const double a[], double b) +{ int j, p, q; + xassert(n >= 2); + /* compute b := b - sum{j : a[j] < 0} */ + for (j = 1; j <= n; j++) + { if (a[j] < 0) + b -= a[j]; + } + /* find a[p] = min{j in 1..n} a[j] */ + p = 1; + for (j = 2; j <= n; j++) + { if (fabs(a[p]) > fabs(a[j])) + p = j; + } + /* find a[q] = min{j in 1..n : j != p} a[j] */ + q = 0; + for (j = 1; j <= n; j++) + { if (j != p) + { if (q == 0 || fabs(a[q]) > fabs(a[j])) + q = j; + } + } + xassert(q != 0); + /* check condition a[p] + a[q] > b */ + return fabs(a[p]) + fabs(a[q]) > b + .001 * (1 + fabs(b)); +} + +/*********************************************************************** +* process_ineq - basic inequality processing +* +* This routine performs basic processing of an inequality of "not +* greater" type +* +* sum{j in 1..n} a[j]*x[j] <= b +* +* specified by the parameters, n, ind, a, and b. +* +* If the inequality can be transformed to "0-1 knapsack" ineqiality +* suitable for generating cover cuts, the routine adds it to the set +* of "0-1 knapsack" inequalities. +* +* Note that the arrays ind and a are not saved on exit. */ + +static void process_ineq(struct csa *csa, int n, int ind[], double a[], + double b, FVS *v) +{ int i; + /* attempt to transform the specified inequality to equivalent or + * relaxed "0-1 knapsack" inequality */ + n = build_ks(csa, n, ind, a, &b, v); + if (n <= 1) + { /* uninteresting inequality (in principle, such inequalities + * should be removed by the preprocessor) */ + goto done; + } + if (!can_be_active(n, a, b)) + { /* inequality is redundant (i.e. cannot be active) */ + goto done; + } + if (is_sos_ineq(n, a, b)) + { /* packing (SOS) inequality is useless for generating cover + * cuts; currently such inequalities are just ignored */ + goto done; + } + /* add resulting "0-1 knapsack" inequality to the set */ + i = glp_add_rows(csa->set, 1); + glp_set_mat_row(csa->set, i, n, ind, a); + glp_set_row_bnds(csa->set, i, GLP_UP, b, b); +done: return; +} + +/**********************************************************************/ + +glp_cov *glp_cov_init(glp_prob *P) +{ /* create and initialize cover cut generator */ + glp_cov *cov; + struct csa csa; + int i, k, len, *ind; + double rhs, *val; + FVS fvs; + csa.P = P; + csa.l = talloc(1+P->n, struct bnd); + csa.u = talloc(1+P->n, struct bnd); + csa.set = glp_create_prob(); + glp_add_cols(csa.set, P->n); + /* initialize bounds of variables with simple bounds */ + init_bounds(&csa); + /* obtain and set variable bounds */ + obtain_vbs(&csa); + /* allocate working arrays */ + ind = talloc(1+P->n, int); + val = talloc(1+P->n, double); + fvs_alloc_vec(&fvs, P->n); + /* process all rows of the root mip */ + for (i = 1; i <= P->m; i++) + { switch (P->row[i]->type) + { case GLP_FR: + break; + case GLP_LO: + /* obtain row of ">=" type */ + len = glp_get_mat_row(P, i, ind, val); + rhs = P->row[i]->lb; + /* transforms it to row of "<=" type */ + for (k = 1; k <= len; k++) + val[k] = - val[k]; + rhs = - rhs; + /* process the row */ + process_ineq(&csa, len, ind, val, rhs, &fvs); + break; + case GLP_UP: + /* obtain row of "<=" type */ + len = glp_get_mat_row(P, i, ind, val); + rhs = P->row[i]->ub; + /* and process it */ + process_ineq(&csa, len, ind, val, rhs, &fvs); + break; + case GLP_DB: + case GLP_FX: + /* double-sided inequalitiy and equality constraints are + * processed as two separate inequalities */ + /* obtain row as if it were of ">=" type */ + len = glp_get_mat_row(P, i, ind, val); + rhs = P->row[i]->lb; + /* transforms it to row of "<=" type */ + for (k = 1; k <= len; k++) + val[k] = - val[k]; + rhs = - rhs; + /* and process it */ + process_ineq(&csa, len, ind, val, rhs, &fvs); + /* obtain the same row as if it were of "<=" type */ + len = glp_get_mat_row(P, i, ind, val); + rhs = P->row[i]->ub; + /* and process it */ + process_ineq(&csa, len, ind, val, rhs, &fvs); + break; + default: + xassert(P != P); + } + } + /* free working arrays */ + tfree(ind); + tfree(val); + fvs_check_vec(&fvs); + fvs_free_vec(&fvs); + /* the set of "0-1 knapsack" inequalities has been built */ + if (csa.set->m == 0) + { /* the set is empty */ + xprintf("No 0-1 knapsack inequalities detected\n"); + cov = NULL; + glp_delete_prob(csa.set); + } + else + { /* create the cover cut generator working area */ + xprintf("Number of 0-1 knapsack inequalities = %d\n", + csa.set->m); + cov = talloc(1, glp_cov); + cov->n = P->n; + cov->set = csa.set; +#if 0 + glp_write_lp(cov->set, 0, "set.lp"); +#endif + } + tfree(csa.l); + tfree(csa.u); + return cov; +} + +/*********************************************************************** +* solve_ks - solve 0-1 knapsack problem +* +* This routine finds (sub)optimal solution to 0-1 knapsack problem: +* +* maximize z = sum{j in 1..n} c[j]x[j] (1) +* +* s.t. sum{j in 1..n} a[j]x[j] <= b (2) +* +* x[j] in {0, 1} for all j in 1..n (3) +* +* It is assumed that the instance is non-normalized, i.e. parameters +* a, b, and c may have any sign. +* +* On exit the routine stores the (sub)optimal point found in locations +* x[1], ..., x[n] and returns the optimal objective value. However, if +* the instance is infeasible, the routine returns INT_MIN. */ + +static int solve_ks(int n, const int a[], int b, const int c[], + char x[]) +{ int z; + /* surprisingly, even for some small instances (n = 50-100) + * MT1 routine takes too much time, so it is used only for tiny + * instances */ + if (n <= 16) +#if 0 + z = ks_enum(n, a, b, c, x); +#else + z = ks_mt1(n, a, b, c, x); +#endif + else + z = ks_greedy(n, a, b, c, x); + return z; +} + +/*********************************************************************** +* simple_cover - find simple cover cut +* +* Given a 0-1 knapsack inequality (which may be globally as well as +* locally valid) +* +* sum{j in 1..n} a[j]x[j] <= b, (1) +* +* where all x[j] are binary variables and all a[j] are positive, and +* a fractional point x~{j in 1..n}, which is feasible to LP relaxation +* of (1), this routine attempts to find a simple cover inequality +* +* sum{j in C} (1 - x[j]) >= 1, (2) +* +* which is valid for (1) and violated at x~. +* +* Actually, the routine finds a cover C, i.e. a subset of {1, ..., n} +* such that +* +* sum{j in C} a[j] > b, (3) +* +* and which minimizes the left-hand side of (2) at x~ +* +* zeta = sum{j in C} (1 - x~[j]). (4) +* +* On exit the routine stores the characteritic vector z{j in 1..n} +* of the cover found (i.e. z[j] = 1 means j in C, and z[j] = 0 means +* j not in C), and returns corresponding minimal value of zeta (4). +* However, if no cover is found, the routine returns DBL_MAX. +* +* ALGORITHM +* +* The separation problem (3)-(4) is converted to 0-1 knapsack problem +* as follows. +* +* First, note that the constraint (3) is equivalent to +* +* sum{j in 1..n} a[j]z[j] >= b + eps, (5) +* +* where eps > 0 is a sufficiently small number (in case of integral +* a and b we may take eps = 1). Multiplying both sides of (5) by (-1) +* gives +* +* sum{j in 1..n} (-a[j])z[j] <= - b - eps. (6) +* +* To make all coefficients in (6) positive, z[j] is complemented by +* substitution z[j] = 1 - z'[j] that finally gives +* +* sum{j in 1..n} a[j]z'[j] <= sum{j in 1..n} a[j] - b - eps. (7) +* +* Minimization of zeta (4) is equivalent to maximization of +* +* -zeta = sum{j in 1..n} (x~[j] - 1)z[j]. (8) +* +* Substitution z[j] = 1 - z'[j] gives +* +* -zeta = sum{j in 1..n} (1 - x~[j])z'[j] - zeta0, (9) +* +* where zeta0 = sum{j in 1..n} (1 - x~[j]) is a constant term. +* +* Thus, the 0-1 knapsack problem to be solved is the following: +* +* maximize +* +* -zeta = sum{j in 1..n} (1 - x~[j])z'[j] - zeta0 (10) +* +* subject to +* +* sum{j in 1..n} a[j]z'[j] <= sum{j in 1..n} a[j] - b - eps (11) +* +* z'[j] in {0,1} for all j = 1,...,n (12) +* +* (The constant term zeta0 doesn't affect the solution, so it can be +* dropped.) */ + +static double simple_cover(int n, const double a[], double b, const + double x[], char z[]) +{ int j, *aa, bb, *cc; + double max_aj, min_aj, s, eps; + xassert(n >= 3); + /* allocate working arrays */ + aa = talloc(1+n, int); + cc = talloc(1+n, int); + /* compute max{j in 1..n} a[j] and min{j in 1..n} a[j] */ + max_aj = 0, min_aj = DBL_MAX; + for (j = 1; j <= n; j++) + { xassert(a[j] > 0); + if (max_aj < a[j]) + max_aj = a[j]; + if (min_aj > a[j]) + min_aj = a[j]; + } + /* scale and round constraint parameters to make them integral; + * note that we make the resulting inequality stronger than (11), + * so a[j]'s are rounded up while rhs is rounded down */ + s = 0; + for (j = 1; j <= n; j++) + { s += a[j]; + aa[j] = ceil(a[j] / max_aj * 1000); + } + bb = floor((s - b) / max_aj * 1000) - 1; + /* scale and round obj. coefficients to make them integral; + * again we make the objective function stronger than (10), so + * the coefficients are rounded down */ + for (j = 1; j <= n; j++) + { xassert(0 <= x[j] && x[j] <= 1); + cc[j] = floor((1 - x[j]) * 1000); + } + /* solve separation problem */ + if (solve_ks(n, aa, bb, cc, z) == INT_MIN) + { /* no cover exists */ + s = DBL_MAX; + goto skip; + } + /* determine z[j] = 1 - z'[j] */ + for (j = 1; j <= n; j++) + { xassert(z[j] == 0 || z[j] == 1); + z[j] ^= 1; + } + /* check condition (11) for original (non-scaled) parameters */ + s = 0; + for (j = 1; j <= n; j++) + { if (z[j]) + s += a[j]; + } + eps = 0.01 * (min_aj >= 1 ? min_aj : 1); + if (!(s >= b + eps)) + { /* no cover found within a precision req'd */ + s = DBL_MAX; + goto skip; + } + /* compute corresponding zeta (4) for cover found */ + s = 0; + for (j = 1; j <= n; j++) + { if (z[j]) + s += 1 - x[j]; + } +skip: /* free working arrays */ + tfree(aa); + tfree(cc); + return s; +} + +/**********************************************************************/ + +void glp_cov_gen1(glp_prob *P, glp_cov *cov, glp_prob *pool) +{ /* generate locally valid simple cover cuts */ + int i, k, len, new_len, *ind; + double *val, rhs, *x, zeta; + char *z; + xassert(P->n == cov->n && P->n == cov->set->n); + xassert(glp_get_status(P) == GLP_OPT); + /* allocate working arrays */ + ind = talloc(1+P->n, int); + val = talloc(1+P->n, double); + x = talloc(1+P->n, double); + z = talloc(1+P->n, char); + /* walk thru 0-1 knapsack inequalities */ + for (i = 1; i <= cov->set->m; i++) + { /* retrieve 0-1 knapsack inequality */ + len = glp_get_mat_row(cov->set, i, ind, val); + rhs = glp_get_row_ub(cov->set, i); + xassert(rhs != +DBL_MAX); + /* FIXME: skip, if slack is too large? */ + /* substitute and eliminate binary variables which have been + * fixed in the current subproblem (this makes the inequality + * only locally valid) */ + new_len = 0; + for (k = 1; k <= len; k++) + { if (glp_get_col_type(P, ind[k]) == GLP_FX) + rhs -= val[k] * glp_get_col_prim(P, ind[k]); + else + { new_len++; + ind[new_len] = ind[k]; + val[new_len] = val[k]; + } + } + len = new_len; + /* we need at least 3 binary variables in the inequality */ + if (len <= 2) + continue; + /* obtain values of binary variables from optimal solution to + * LP relaxation of current subproblem */ + for (k = 1; k <= len; k++) + { xassert(glp_get_col_kind(P, ind[k]) == GLP_BV); + x[k] = glp_get_col_prim(P, ind[k]); + if (x[k] < 0.00001) + x[k] = 0; + else if (x[k] > 0.99999) + x[k] = 1; + /* if val[k] < 0, perform substitution x[k] = 1 - x'[k] to + * make all coefficients positive */ + if (val[k] < 0) + { ind[k] = - ind[k]; /* x[k] is complemented */ + val[k] = - val[k]; + rhs += val[k]; + x[k] = 1 - x[k]; + } + } + /* find locally valid simple cover cut */ + zeta = simple_cover(len, val, rhs, x, z); + if (zeta > 0.95) + { /* no violation or insufficient violation; see (2) */ + continue; + } + /* construct cover inequality (2) for the cover found, which + * for original binary variables x[k] is equivalent to: + * sum{k in C'} x[k] + sum{k in C"} x'[k] <= |C| - 1 + * or + * sum{k in C'} x[k] + sum{k in C"} (1 - x[k]) <= |C| - 1 + * or + * sum{k in C'} x[k] - sum{k in C"} x[k] <= |C'| - 1 + * since |C| - |C"| = |C'| */ + new_len = 0; + rhs = -1; + for (k = 1; k <= len; k++) + { if (z[k]) + { new_len++; + if (ind[k] > 0) + { ind[new_len] = +ind[k]; + val[new_len] = +1; + rhs++; + } + else /* ind[k] < 0 */ + { ind[new_len] = -ind[k]; + val[new_len] = -1; + } + } + } + len = new_len; + /* add the cover inequality to the local cut pool */ + k = glp_add_rows(pool, 1); + glp_set_mat_row(pool, k, len, ind, val); + glp_set_row_bnds(pool, k, GLP_UP, rhs, rhs); + } + /* free working arrays */ + tfree(ind); + tfree(val); + tfree(x); + tfree(z); + return; +} + +/**********************************************************************/ + +void glp_cov_free(glp_cov *cov) +{ /* delete cover cut generator workspace */ + xassert(cov != NULL); + glp_delete_prob(cov->set); + tfree(cov); + return; +} + +/* eof */ diff --git a/test/monniaux/glpk-4.65/src/intopt/fpump.c b/test/monniaux/glpk-4.65/src/intopt/fpump.c new file mode 100644 index 00000000..0bdd6d4e --- /dev/null +++ b/test/monniaux/glpk-4.65/src/intopt/fpump.c @@ -0,0 +1,360 @@ +/* fpump.c (feasibility pump heuristic) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2009-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 "env.h" +#include "ios.h" +#include "rng.h" + +/*********************************************************************** +* NAME +* +* ios_feas_pump - feasibility pump heuristic +* +* SYNOPSIS +* +* #include "glpios.h" +* void ios_feas_pump(glp_tree *T); +* +* DESCRIPTION +* +* The routine ios_feas_pump is a simple implementation of the Feasi- +* bility Pump heuristic. +* +* REFERENCES +* +* M.Fischetti, F.Glover, and A.Lodi. "The feasibility pump." Math. +* Program., Ser. A 104, pp. 91-104 (2005). */ + +struct VAR +{ /* binary variable */ + int j; + /* ordinal number */ + int x; + /* value in the rounded solution (0 or 1) */ + double d; + /* sorting key */ +}; + +static int CDECL fcmp(const void *x, const void *y) +{ /* comparison routine */ + const struct VAR *vx = x, *vy = y; + if (vx->d > vy->d) + return -1; + else if (vx->d < vy->d) + return +1; + else + return 0; +} + +void ios_feas_pump(glp_tree *T) +{ glp_prob *P = T->mip; + int n = P->n; + glp_prob *lp = NULL; + struct VAR *var = NULL; + RNG *rand = NULL; + GLPCOL *col; + glp_smcp parm; + int j, k, new_x, nfail, npass, nv, ret, stalling; + double dist, tol; + xassert(glp_get_status(P) == GLP_OPT); + /* this heuristic is applied only once on the root level */ + if (!(T->curr->level == 0 && T->curr->solved == 1)) goto done; + /* determine number of binary variables */ + nv = 0; + for (j = 1; j <= n; j++) + { col = P->col[j]; + /* if x[j] is continuous, skip it */ + if (col->kind == GLP_CV) continue; + /* if x[j] is fixed, skip it */ + if (col->type == GLP_FX) continue; + /* x[j] is non-fixed integer */ + xassert(col->kind == GLP_IV); + if (col->type == GLP_DB && col->lb == 0.0 && col->ub == 1.0) + { /* x[j] is binary */ + nv++; + } + else + { /* x[j] is general integer */ + if (T->parm->msg_lev >= GLP_MSG_ALL) + xprintf("FPUMP heuristic cannot be applied due to genera" + "l integer variables\n"); + goto done; + } + } + /* there must be at least one binary variable */ + if (nv == 0) goto done; + if (T->parm->msg_lev >= GLP_MSG_ALL) + xprintf("Applying FPUMP heuristic...\n"); + /* build the list of binary variables */ + var = xcalloc(1+nv, sizeof(struct VAR)); + k = 0; + for (j = 1; j <= n; j++) + { col = P->col[j]; + if (col->kind == GLP_IV && col->type == GLP_DB) + var[++k].j = j; + } + xassert(k == nv); + /* create working problem object */ + lp = glp_create_prob(); +more: /* copy the original problem object to keep it intact */ + glp_copy_prob(lp, P, GLP_OFF); + /* we are interested to find an integer feasible solution, which + is better than the best known one */ + if (P->mip_stat == GLP_FEAS) + { int *ind; + double *val, bnd; + /* add a row and make it identical to the objective row */ + glp_add_rows(lp, 1); + ind = xcalloc(1+n, sizeof(int)); + val = xcalloc(1+n, sizeof(double)); + for (j = 1; j <= n; j++) + { ind[j] = j; + val[j] = P->col[j]->coef; + } + glp_set_mat_row(lp, lp->m, n, ind, val); + xfree(ind); + xfree(val); + /* introduce upper (minimization) or lower (maximization) + bound to the original objective function; note that this + additional constraint is not violated at the optimal point + to LP relaxation */ +#if 0 /* modified by xypron <xypron.glpk@gmx.de> */ + if (P->dir == GLP_MIN) + { bnd = P->mip_obj - 0.10 * (1.0 + fabs(P->mip_obj)); + if (bnd < P->obj_val) bnd = P->obj_val; + glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0); + } + else if (P->dir == GLP_MAX) + { bnd = P->mip_obj + 0.10 * (1.0 + fabs(P->mip_obj)); + if (bnd > P->obj_val) bnd = P->obj_val; + glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0); + } + else + xassert(P != P); +#else + bnd = 0.1 * P->obj_val + 0.9 * P->mip_obj; + /* xprintf("bnd = %f\n", bnd); */ + if (P->dir == GLP_MIN) + glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0); + else if (P->dir == GLP_MAX) + glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0); + else + xassert(P != P); +#endif + } + /* reset pass count */ + npass = 0; + /* invalidate the rounded point */ + for (k = 1; k <= nv; k++) + var[k].x = -1; +pass: /* next pass starts here */ + npass++; + if (T->parm->msg_lev >= GLP_MSG_ALL) + xprintf("Pass %d\n", npass); + /* initialize minimal distance between the basic point and the + rounded one obtained during this pass */ + dist = DBL_MAX; + /* reset failure count (the number of succeeded iterations failed + to improve the distance) */ + nfail = 0; + /* if it is not the first pass, perturb the last rounded point + rather than construct it from the basic solution */ + if (npass > 1) + { double rho, temp; + if (rand == NULL) + rand = rng_create_rand(); + for (k = 1; k <= nv; k++) + { j = var[k].j; + col = lp->col[j]; + rho = rng_uniform(rand, -0.3, 0.7); + if (rho < 0.0) rho = 0.0; + temp = fabs((double)var[k].x - col->prim); + if (temp + rho > 0.5) var[k].x = 1 - var[k].x; + } + goto skip; + } +loop: /* innermost loop begins here */ + /* round basic solution (which is assumed primal feasible) */ + stalling = 1; + for (k = 1; k <= nv; k++) + { col = lp->col[var[k].j]; + if (col->prim < 0.5) + { /* rounded value is 0 */ + new_x = 0; + } + else + { /* rounded value is 1 */ + new_x = 1; + } + if (var[k].x != new_x) + { stalling = 0; + var[k].x = new_x; + } + } + /* if the rounded point has not changed (stalling), choose and + flip some its entries heuristically */ + if (stalling) + { /* compute d[j] = |x[j] - round(x[j])| */ + for (k = 1; k <= nv; k++) + { col = lp->col[var[k].j]; + var[k].d = fabs(col->prim - (double)var[k].x); + } + /* sort the list of binary variables by descending d[j] */ + qsort(&var[1], nv, sizeof(struct VAR), fcmp); + /* choose and flip some rounded components */ + for (k = 1; k <= nv; k++) + { if (k >= 5 && var[k].d < 0.35 || k >= 10) break; + var[k].x = 1 - var[k].x; + } + } +skip: /* check if the time limit has been exhausted */ + if (T->parm->tm_lim < INT_MAX && + (double)(T->parm->tm_lim - 1) <= + 1000.0 * xdifftime(xtime(), T->tm_beg)) goto done; + /* build the objective, which is the distance between the current + (basic) point and the rounded one */ + lp->dir = GLP_MIN; + lp->c0 = 0.0; + for (j = 1; j <= n; j++) + lp->col[j]->coef = 0.0; + for (k = 1; k <= nv; k++) + { j = var[k].j; + if (var[k].x == 0) + lp->col[j]->coef = +1.0; + else + { lp->col[j]->coef = -1.0; + lp->c0 += 1.0; + } + } + /* minimize the distance with the simplex method */ + glp_init_smcp(&parm); + if (T->parm->msg_lev <= GLP_MSG_ERR) + parm.msg_lev = T->parm->msg_lev; + else if (T->parm->msg_lev <= GLP_MSG_ALL) + { parm.msg_lev = GLP_MSG_ON; + parm.out_dly = 10000; + } + ret = glp_simplex(lp, &parm); + if (ret != 0) + { if (T->parm->msg_lev >= GLP_MSG_ERR) + xprintf("Warning: glp_simplex returned %d\n", ret); + goto done; + } + ret = glp_get_status(lp); + if (ret != GLP_OPT) + { if (T->parm->msg_lev >= GLP_MSG_ERR) + xprintf("Warning: glp_get_status returned %d\n", ret); + goto done; + } + if (T->parm->msg_lev >= GLP_MSG_DBG) + xprintf("delta = %g\n", lp->obj_val); + /* check if the basic solution is integer feasible; note that it + may be so even if the minimial distance is positive */ + tol = 0.3 * T->parm->tol_int; + for (k = 1; k <= nv; k++) + { col = lp->col[var[k].j]; + if (tol < col->prim && col->prim < 1.0 - tol) break; + } + if (k > nv) + { /* okay; the basic solution seems to be integer feasible */ + double *x = xcalloc(1+n, sizeof(double)); + for (j = 1; j <= n; j++) + { x[j] = lp->col[j]->prim; + if (P->col[j]->kind == GLP_IV) x[j] = floor(x[j] + 0.5); + } +#if 1 /* modified by xypron <xypron.glpk@gmx.de> */ + /* reset direction and right-hand side of objective */ + lp->c0 = P->c0; + lp->dir = P->dir; + /* fix integer variables */ + for (k = 1; k <= nv; k++) +#if 0 /* 18/VI-2013; fixed by mao + * this bug causes numerical instability, because column statuses + * are not changed appropriately */ + { lp->col[var[k].j]->lb = x[var[k].j]; + lp->col[var[k].j]->ub = x[var[k].j]; + lp->col[var[k].j]->type = GLP_FX; + } +#else + glp_set_col_bnds(lp, var[k].j, GLP_FX, x[var[k].j], 0.); +#endif + /* copy original objective function */ + for (j = 1; j <= n; j++) + lp->col[j]->coef = P->col[j]->coef; + /* solve original LP and copy result */ + ret = glp_simplex(lp, &parm); + if (ret != 0) + { if (T->parm->msg_lev >= GLP_MSG_ERR) + xprintf("Warning: glp_simplex returned %d\n", ret); +#if 1 /* 17/III-2016: fix memory leak */ + xfree(x); +#endif + goto done; + } + ret = glp_get_status(lp); + if (ret != GLP_OPT) + { if (T->parm->msg_lev >= GLP_MSG_ERR) + xprintf("Warning: glp_get_status returned %d\n", ret); +#if 1 /* 17/III-2016: fix memory leak */ + xfree(x); +#endif + goto done; + } + for (j = 1; j <= n; j++) + if (P->col[j]->kind != GLP_IV) x[j] = lp->col[j]->prim; +#endif + ret = glp_ios_heur_sol(T, x); + xfree(x); + if (ret == 0) + { /* the integer solution is accepted */ + if (ios_is_hopeful(T, T->curr->bound)) + { /* it is reasonable to apply the heuristic once again */ + goto more; + } + else + { /* the best known integer feasible solution just found + is close to optimal solution to LP relaxation */ + goto done; + } + } + } + /* the basic solution is fractional */ + if (dist == DBL_MAX || + lp->obj_val <= dist - 1e-6 * (1.0 + dist)) + { /* the distance is reducing */ + nfail = 0, dist = lp->obj_val; + } + else + { /* improving the distance failed */ + nfail++; + } + if (nfail < 3) goto loop; + if (npass < 5) goto pass; +done: /* delete working objects */ + if (lp != NULL) glp_delete_prob(lp); + if (var != NULL) xfree(var); + if (rand != NULL) rng_delete_rand(rand); + return; +} + +/* eof */ diff --git a/test/monniaux/glpk-4.65/src/intopt/gmicut.c b/test/monniaux/glpk-4.65/src/intopt/gmicut.c new file mode 100644 index 00000000..4ef0b746 --- /dev/null +++ b/test/monniaux/glpk-4.65/src/intopt/gmicut.c @@ -0,0 +1,284 @@ +/* gmicut.c (Gomory's mixed integer cut generator) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2002-2016 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 "env.h" +#include "prob.h" + +/*********************************************************************** +* NAME +* +* glp_gmi_cut - generate Gomory's mixed integer cut (core routine) +* +* SYNOPSIS +* +* int glp_gmi_cut(glp_prob *P, int j, int ind[], double val[], double +* phi[]); +* +* DESCRIPTION +* +* This routine attempts to generate a Gomory's mixed integer cut for +* specified integer column (structural variable), whose primal value +* in current basic solution is integer infeasible (fractional). +* +* On entry to the routine the basic solution contained in the problem +* object P should be optimal, and the basis factorization should be +* valid. The parameter j should specify the ordinal number of column +* (structural variable x[j]), for which the cut should be generated, +* 1 <= j <= n, where n is the number of columns in the problem object. +* This column should be integer, non-fixed, and basic, and its primal +* value should be fractional. +* +* The cut generated by the routine is the following inequality: +* +* sum a[j] * x[j] >= b, +* +* which is expected to be violated at the current basic solution. +* +* If the cut has been successfully generated, the routine stores its +* non-zero coefficients a[j] and corresponding column indices j in the +* array locations val[1], ..., val[len] and ind[1], ..., ind[len], +* where 1 <= len <= n is the number of non-zero coefficients. The +* right-hand side value b is stored in val[0], and ind[0] is set to 0. +* +* The working array phi should have 1+m+n locations (location phi[0] +* is not used), where m and n is the number of rows and columns in the +* problem object, resp. +* +* RETURNS +* +* If the cut has been successfully generated, the routine returns +* len, the number of non-zero coefficients in the cut, 1 <= len <= n. +* +* Otherwise, the routine returns one of the following codes: +* +* -1 current basis factorization is not valid; +* +* -2 current basic solution is not optimal; +* +* -3 column ordinal number j is out of range; +* +* -4 variable x[j] is not of integral kind; +* +* -5 variable x[j] is either fixed or non-basic; +* +* -6 primal value of variable x[j] in basic solution is too close +* to nearest integer; +* +* -7 some coefficients in the simplex table row corresponding to +* variable x[j] are too large in magnitude; +* +* -8 some free (unbounded) variables have non-zero coefficients in +* the simplex table row corresponding to variable x[j]. +* +* ALGORITHM +* +* See glpk/doc/notes/gomory (in Russian). */ + +#define f(x) ((x) - floor(x)) +/* compute fractional part of x */ + +int glp_gmi_cut(glp_prob *P, int j, + int ind[/*1+n*/], double val[/*1+n*/], double phi[/*1+m+n*/]) +{ int m = P->m; + int n = P->n; + GLPROW *row; + GLPCOL *col; + GLPAIJ *aij; + int i, k, len, kind, stat; + double lb, ub, alfa, beta, ksi, phi1, rhs; + /* sanity checks */ + if (!(P->m == 0 || P->valid)) + { /* current basis factorization is not valid */ + return -1; + } + if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS)) + { /* current basic solution is not optimal */ + return -2; + } + if (!(1 <= j && j <= n)) + { /* column ordinal number is out of range */ + return -3; + } + col = P->col[j]; + if (col->kind != GLP_IV) + { /* x[j] is not of integral kind */ + return -4; + } + if (col->type == GLP_FX || col->stat != GLP_BS) + { /* x[j] is either fixed or non-basic */ + return -5; + } + if (fabs(col->prim - floor(col->prim + 0.5)) < 0.001) + { /* primal value of x[j] is too close to nearest integer */ + return -6; + } + /* compute row of the simplex tableau, which (row) corresponds + * to specified basic variable xB[i] = x[j]; see (23) */ + len = glp_eval_tab_row(P, m+j, ind, val); + /* determine beta[i], which a value of xB[i] in optimal solution + * to current LP relaxation; note that this value is the same as + * if it would be computed with formula (27); it is assumed that + * beta[i] is fractional enough */ + beta = P->col[j]->prim; + /* compute cut coefficients phi and right-hand side rho, which + * correspond to formula (30); dense format is used, because rows + * of the simplex tableau are usually dense */ + for (k = 1; k <= m+n; k++) + phi[k] = 0.0; + rhs = f(beta); /* initial value of rho; see (28), (32) */ + for (j = 1; j <= len; j++) + { /* determine original number of non-basic variable xN[j] */ + k = ind[j]; + xassert(1 <= k && k <= m+n); + /* determine the kind, bounds and current status of xN[j] in + * optimal solution to LP relaxation */ + if (k <= m) + { /* auxiliary variable */ + row = P->row[k]; + kind = GLP_CV; + lb = row->lb; + ub = row->ub; + stat = row->stat; + } + else + { /* structural variable */ + col = P->col[k-m]; + kind = col->kind; + lb = col->lb; + ub = col->ub; + stat = col->stat; + } + /* xN[j] cannot be basic */ + xassert(stat != GLP_BS); + /* determine row coefficient ksi[i,j] at xN[j]; see (23) */ + ksi = val[j]; + /* if ksi[i,j] is too large in magnitude, report failure */ + if (fabs(ksi) > 1e+05) + return -7; + /* if ksi[i,j] is too small in magnitude, skip it */ + if (fabs(ksi) < 1e-10) + goto skip; + /* compute row coefficient alfa[i,j] at y[j]; see (26) */ + switch (stat) + { case GLP_NF: + /* xN[j] is free (unbounded) having non-zero ksi[i,j]; + * report failure */ + return -8; + case GLP_NL: + /* xN[j] has active lower bound */ + alfa = - ksi; + break; + case GLP_NU: + /* xN[j] has active upper bound */ + alfa = + ksi; + break; + case GLP_NS: + /* xN[j] is fixed; skip it */ + goto skip; + default: + xassert(stat != stat); + } + /* compute cut coefficient phi'[j] at y[j]; see (21), (28) */ + switch (kind) + { case GLP_IV: + /* y[j] is integer */ + if (fabs(alfa - floor(alfa + 0.5)) < 1e-10) + { /* alfa[i,j] is close to nearest integer; skip it */ + goto skip; + } + else if (f(alfa) <= f(beta)) + phi1 = f(alfa); + else + phi1 = (f(beta) / (1.0 - f(beta))) * (1.0 - f(alfa)); + break; + case GLP_CV: + /* y[j] is continuous */ + if (alfa >= 0.0) + phi1 = + alfa; + else + phi1 = (f(beta) / (1.0 - f(beta))) * (- alfa); + break; + default: + xassert(kind != kind); + } + /* compute cut coefficient phi[j] at xN[j] and update right- + * hand side rho; see (31), (32) */ + switch (stat) + { case GLP_NL: + /* xN[j] has active lower bound */ + phi[k] = + phi1; + rhs += phi1 * lb; + break; + case GLP_NU: + /* xN[j] has active upper bound */ + phi[k] = - phi1; + rhs -= phi1 * ub; + break; + default: + xassert(stat != stat); + } +skip: ; + } + /* now the cut has the form sum_k phi[k] * x[k] >= rho, where cut + * coefficients are stored in the array phi in dense format; + * x[1,...,m] are auxiliary variables, x[m+1,...,m+n] are struc- + * tural variables; see (30) */ + /* eliminate auxiliary variables in order to express the cut only + * through structural variables; see (33) */ + for (i = 1; i <= m; i++) + { if (fabs(phi[i]) < 1e-10) + continue; + /* auxiliary variable x[i] has non-zero cut coefficient */ + row = P->row[i]; + /* x[i] cannot be fixed variable */ + xassert(row->type != GLP_FX); + /* substitute x[i] = sum_j a[i,j] * x[m+j] */ + for (aij = row->ptr; aij != NULL; aij = aij->r_next) + phi[m+aij->col->j] += phi[i] * aij->val; + } + /* convert the final cut to sparse format and substitute fixed + * (structural) variables */ + len = 0; + for (j = 1; j <= n; j++) + { if (fabs(phi[m+j]) < 1e-10) + continue; + /* structural variable x[m+j] has non-zero cut coefficient */ + col = P->col[j]; + if (col->type == GLP_FX) + { /* eliminate x[m+j] */ + rhs -= phi[m+j] * col->lb; + } + else + { len++; + ind[len] = j; + val[len] = phi[m+j]; + } + } + if (fabs(rhs) < 1e-12) + rhs = 0.0; + ind[0] = 0, val[0] = rhs; + /* the cut has been successfully generated */ + return len; +} + +/* eof */ diff --git a/test/monniaux/glpk-4.65/src/intopt/gmigen.c b/test/monniaux/glpk-4.65/src/intopt/gmigen.c new file mode 100644 index 00000000..627682cb --- /dev/null +++ b/test/monniaux/glpk-4.65/src/intopt/gmigen.c @@ -0,0 +1,142 @@ +/* gmigen.c (Gomory's mixed integer cuts generator) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2002-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 "env.h" +#include "prob.h" + +/*********************************************************************** +* NAME +* +* glp_gmi_gen - generate Gomory's mixed integer cuts +* +* SYNOPSIS +* +* int glp_gmi_gen(glp_prob *P, glp_prob *pool, int max_cuts); +* +* DESCRIPTION +* +* This routine attempts to generate Gomory's mixed integer cuts for +* integer variables, whose primal values in current basic solution are +* integer infeasible (fractional). +* +* On entry to the routine the basic solution contained in the problem +* object P should be optimal, and the basis factorization should be +* valid. +* +* The cutting plane inequalities generated by the routine are added to +* the specified cut pool. +* +* The parameter max_cuts specifies the maximal number of cuts to be +* generated. Note that the number of cuts cannot exceed the number of +* basic variables, which is the number of rows in the problem object. +* +* RETURNS +* +* The routine returns the number of cuts that have been generated and +* added to the cut pool. */ + +#define f(x) ((x) - floor(x)) +/* compute fractional part of x */ + +struct var { int j; double f; }; + +static int CDECL fcmp(const void *p1, const void *p2) +{ const struct var *v1 = p1, *v2 = p2; + if (v1->f > v2->f) return -1; + if (v1->f < v2->f) return +1; + return 0; +} + +int glp_gmi_gen(glp_prob *P, glp_prob *pool, int max_cuts) +{ int m = P->m; + int n = P->n; + GLPCOL *col; + struct var *var; + int i, j, k, t, len, nv, nnn, *ind; + double frac, *val, *phi; + /* sanity checks */ + if (!(P->m == 0 || P->valid)) + xerror("glp_gmi_gen: basis factorization does not exist\n"); + if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS)) + xerror("glp_gmi_gen: optimal basic solution required\n"); + if (pool->n != n) + xerror("glp_gmi_gen: cut pool has wrong number of columns\n"); + /* allocate working arrays */ + var = xcalloc(1+n, sizeof(struct var)); + ind = xcalloc(1+n, sizeof(int)); + val = xcalloc(1+n, sizeof(double)); + phi = xcalloc(1+m+n, sizeof(double)); + /* build the list of integer structural variables, which are + * basic and have integer infeasible (fractional) primal values + * in optimal solution to specified LP */ + nv = 0; + for (j = 1; j <= n; j++) + { col = P->col[j]; + if (col->kind != GLP_IV) + continue; + if (col->type == GLP_FX) + continue; + if (col->stat != GLP_BS) + continue; + frac = f(col->prim); + if (!(0.05 <= frac && frac <= 0.95)) + continue; + /* add variable to the list */ + nv++, var[nv].j = j, var[nv].f = frac; + } + /* sort the list by descending fractionality */ + qsort(&var[1], nv, sizeof(struct var), fcmp); + /* try to generate cuts by one for each variable in the list, but + * not more than max_cuts cuts */ + nnn = 0; + for (t = 1; t <= nv; t++) + { len = glp_gmi_cut(P, var[t].j, ind, val, phi); + if (len < 1) + goto skip; + /* if the cut inequality seems to be badly scaled, reject it + * to avoid numerical difficulties */ + for (k = 1; k <= len; k++) + { if (fabs(val[k]) < 1e-03) + goto skip; + if (fabs(val[k]) > 1e+03) + goto skip; + } + /* add the cut to the cut pool for further consideration */ + i = glp_add_rows(pool, 1); + glp_set_row_bnds(pool, i, GLP_LO, val[0], 0); + glp_set_mat_row(pool, i, len, ind, val); + /* one cut has been generated */ + nnn++; + if (nnn == max_cuts) + break; +skip: ; + } + /* free working arrays */ + xfree(var); + xfree(ind); + xfree(val); + xfree(phi); + return nnn; +} + +/* eof */ diff --git a/test/monniaux/glpk-4.65/src/intopt/mirgen.c b/test/monniaux/glpk-4.65/src/intopt/mirgen.c new file mode 100644 index 00000000..45a0a55d --- /dev/null +++ b/test/monniaux/glpk-4.65/src/intopt/mirgen.c @@ -0,0 +1,1529 @@ +/* mirgen.c (mixed integer rounding cuts generator) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2007-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/>. +***********************************************************************/ + +#if 1 /* 29/II-2016 by Chris */ +/*---------------------------------------------------------------------- +Subject: Mir cut generation performance improvement +From: Chris Matrakidis <cmatraki@gmail.com> +To: Andrew Makhorin <mao@gnu.org>, help-glpk <help-glpk@gnu.org> + +Andrew, + +I noticed that mir cut generation takes considerable time on some large +problems (like rocII-4-11 from miplib). The attached patch makes two +improvements that considerably improve performance in such instances: +1. A lot of time was spent on generating a temporary vector in function +aggregate_row. It is a lot faster to reuse an existing vector. +2. A search for an element in the same function was done in row order, +where using the elements in the order they are in the column is more +efficient. This changes the generated cuts in some cases, but seems +neutral overall (0.3% less cuts in a test set of 64 miplib instances). + +Best Regards, + +Chris Matrakidis +----------------------------------------------------------------------*/ +#endif + +#include "env.h" +#include "prob.h" +#include "spv.h" + +#define MIR_DEBUG 0 + +#define MAXAGGR 5 +/* maximal number of rows that can be aggregated */ + +struct glp_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 */ + SPV *agg_vec; /* SPV 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 */ + SPV *mod_vec; /* SPV 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 */ + SPV *cut_vec; /* SPV 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 +* +* glp_mir_init - create and initialize MIR cut generator +* +* SYNOPSIS +* +* glp_mir *glp_mir_init(glp_prob *P); +* +* DESCRIPTION +* +* This routine creates and initializes the MIR cut generator for the +* specified problem object. +* +* RETURNS +* +* The routine returns a pointer to the MIR cut generator workspace. */ + +static void set_row_attrib(glp_prob *mip, glp_mir *mir) +{ /* set global row attributes */ + 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_prob *mip, glp_mir *mir) +{ /* set global column attributes */ + 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_prob *mip, glp_mir *mir) +{ /* set variable bounds */ + 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_prob *mip, glp_mir *mir) +{ /* mark rows which should not be used */ + 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; +} + +glp_mir *glp_mir_init(glp_prob *mip) +{ /* create and initialize MIR cut generator */ + int m = mip->m; + int n = mip->n; + glp_mir *mir; +#if MIR_DEBUG + xprintf("ios_mir_init: warning: debug mode enabled\n"); +#endif + /* allocate working area */ + mir = xmalloc(sizeof(glp_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 = spv_create_vec(m+n); + mir->subst = xcalloc(1+m+n, sizeof(char)); + mir->mod_vec = spv_create_vec(m+n); + mir->cut_vec = spv_create_vec(m+n); + /* set global row attributes */ + set_row_attrib(mip, mir); + /* set global column attributes */ + set_col_attrib(mip, mir); + /* set variable bounds */ + set_var_bounds(mip, mir); + /* mark rows which should not be used */ + mark_useless_rows(mip, mir); + return mir; +} + +/*********************************************************************** +* NAME +* +* glp_mir_gen - generate mixed integer rounding (MIR) cuts +* +* SYNOPSIS +* +* int glp_mir_gen(glp_prob *P, glp_mir *mir, glp_prob *pool); +* +* DESCRIPTION +* +* This routine attempts to generate mixed integer rounding (MIR) cuts +* for current basic solution to the specified problem object. +* +* The cutting plane inequalities generated by the routine are added to +* the specified cut pool. +* +* RETURNS +* +* The routine returns the number of cuts that have been generated and +* added to the cut pool. */ + +static void get_current_point(glp_prob *mip, glp_mir *mir) +{ /* obtain current point */ + 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(glp_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_prob *mip, glp_mir *mir, int i) +{ /* use original i-th row as initial aggregated constraint */ + 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 */ + spv_clear_vec(mir->agg_vec); + spv_set_vj(mir->agg_vec, i, 1.0); + for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next) + spv_set_vj(mir->agg_vec, m + aij->col->j, - aij->val); + mir->agg_rhs = 0.0; +#if MIR_DEBUG + spv_check_vec(mir->agg_vec); +#endif + return; +} + +#if MIR_DEBUG +static void check_agg_row(glp_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(glp_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 */ + spv_clean_vec(mir->agg_vec, DBL_EPSILON); +#if MIR_DEBUG + spv_check_vec(mir->agg_vec); +#endif + return; +} + +static void bound_subst_heur(glp_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(glp_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 */ + spv_copy_vec(mir->mod_vec, mir->agg_vec); + mir->mod_rhs = mir->agg_rhs; +#if MIR_DEBUG + spv_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) + { spv_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) + { spv_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 + spv_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 + spv_check_vec(mir->mod_vec); +#endif + return; +} + +#if MIR_DEBUG +static void check_mod_row(glp_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 CDECL 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(glp_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; + spv_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 */ + spv_clean_vec(mir->cut_vec, DBL_EPSILON); +#if MIR_DEBUG + spv_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; + } + spv_clean_vec(mir->cut_vec, 0.0); +#if MIR_DEBUG + spv_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 + spv_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); +#if 0 /* 06/III-2016; notorious bug reported many times */ + xassert(x[j] >= -0.001); +#else + if (x[j] < -0.001) + { xprintf("glp_mir_gen: warning: x[%d] = %g\n", j, x[j]); + r_best = 0.0; + goto skip; + } +#endif + 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); +#if 0 /* 06/III-2016; notorious bug reported many times */ + xassert(x >= -0.001); +#else + if (x < -0.001) + { xprintf("glp_mir_gen: warning: x = %g\n", x); + r_best = 0.0; + goto skip; + } +#endif + 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 + spv_check_vec(mir->cut_vec); +#endif +skip: /* free working arrays */ + xfree(u); + xfree(x); + xfree(alpha); +done: return r_best; +} + +#if MIR_DEBUG +static void check_raw_cut(glp_mir *mir, double r_best) +{ /* check raw cut before back bound substitution */ + 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->cut_vec->nnz; j++) + { k = mir->cut_vec->ind[j]; + xassert(1 <= k && k <= m+n); + 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); + r += mir->cut_vec->val[j] * x; + if (big < fabs(mir->cut_vec->val[j])) + big = fabs(mir->cut_vec->val[j]); + } + r -= mir->cut_rhs; + if (big < fabs(mir->cut_rhs)) + big = fabs(mir->cut_rhs); + /* the residual must be close to r_best */ + xassert(fabs(r - r_best) <= 1e-6 * big); + return; +} +#endif + +static void back_subst(glp_mir *mir) +{ /* back substitution of original bounds */ + int m = mir->m; + int n = mir->n; + int j, jj, k, kk; + /* at first, restore bounds of integer variables (because on + restoring variable bounds of continuous variables we need + original, not shifted, bounds of integer variables) */ + 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]) continue; /* skip continuous */ + if (mir->subst[k] == 'L') + { /* x'[k] = x[k] - lb[k] */ + xassert(mir->lb[k] != -DBL_MAX); + xassert(mir->vlb[k] == 0); + mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k]; + } + else if (mir->subst[k] == 'U') + { /* x'[k] = ub[k] - x[k] */ + xassert(mir->ub[k] != +DBL_MAX); + xassert(mir->vub[k] == 0); + mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k]; + mir->cut_vec->val[j] = - mir->cut_vec->val[j]; + } + else + xassert(k != k); + } + /* now restore bounds of continuous variables */ + 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]) continue; /* skip integer */ + 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'[k] = x[k] - lb[k] */ + mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k]; + } + else + { /* x'[k] = x[k] - lb[k] * x[kk] */ + jj = mir->cut_vec->pos[kk]; +#if 0 + xassert(jj != 0); +#else + if (jj == 0) + { spv_set_vj(mir->cut_vec, kk, 1.0); + jj = mir->cut_vec->pos[kk]; + xassert(jj != 0); + mir->cut_vec->val[jj] = 0.0; + } +#endif + mir->cut_vec->val[jj] -= mir->cut_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->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k]; + } + else + { /* x'[k] = ub[k] * x[kk] - x[k] */ + jj = mir->cut_vec->pos[kk]; + if (jj == 0) + { spv_set_vj(mir->cut_vec, kk, 1.0); + jj = mir->cut_vec->pos[kk]; + xassert(jj != 0); + mir->cut_vec->val[jj] = 0.0; + } + mir->cut_vec->val[jj] += mir->cut_vec->val[j] * + mir->ub[k]; + } + mir->cut_vec->val[j] = - mir->cut_vec->val[j]; + } + else + xassert(k != k); + } +#if MIR_DEBUG + spv_check_vec(mir->cut_vec); +#endif + return; +} + +#if MIR_DEBUG +static void check_cut_row(glp_mir *mir, double r_best) +{ /* check the cut after back bound substitution or elimination of + auxiliary variables */ + 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->cut_vec->nnz; j++) + { k = mir->cut_vec->ind[j]; + xassert(1 <= k && k <= m+n); + r += mir->cut_vec->val[j] * mir->x[k]; + if (big < fabs(mir->cut_vec->val[j])) + big = fabs(mir->cut_vec->val[j]); + } + r -= mir->cut_rhs; + if (big < fabs(mir->cut_rhs)) + big = fabs(mir->cut_rhs); + /* the residual must be close to r_best */ + xassert(fabs(r - r_best) <= 1e-6 * big); + return; +} +#endif + +static void subst_aux_vars(glp_prob *mip, glp_mir *mir) +{ /* final substitution to eliminate auxiliary variables */ + int m = mir->m; + int n = mir->n; + GLPAIJ *aij; + int j, k, kk, jj; + for (j = mir->cut_vec->nnz; j >= 1; j--) + { k = mir->cut_vec->ind[j]; + xassert(1 <= k && k <= m+n); + if (k > m) continue; /* skip structurals */ + for (aij = mip->row[k]->ptr; aij != NULL; aij = aij->r_next) + { kk = m + aij->col->j; /* structural */ + jj = mir->cut_vec->pos[kk]; + if (jj == 0) + { spv_set_vj(mir->cut_vec, kk, 1.0); + jj = mir->cut_vec->pos[kk]; + mir->cut_vec->val[jj] = 0.0; + } + mir->cut_vec->val[jj] += mir->cut_vec->val[j] * aij->val; + } + mir->cut_vec->val[j] = 0.0; + } + spv_clean_vec(mir->cut_vec, 0.0); + return; +} + +static void add_cut(glp_mir *mir, glp_prob *pool) +{ /* add constructed cut inequality to the cut pool */ + int m = mir->m; + int n = mir->n; + int j, k, len; + int *ind = xcalloc(1+n, sizeof(int)); + double *val = xcalloc(1+n, sizeof(double)); + len = 0; + for (j = mir->cut_vec->nnz; j >= 1; j--) + { k = mir->cut_vec->ind[j]; + xassert(m+1 <= k && k <= m+n); + len++, ind[len] = k - m, val[len] = mir->cut_vec->val[j]; + } +#if 0 +#if 0 + ios_add_cut_row(tree, pool, GLP_RF_MIR, len, ind, val, GLP_UP, + mir->cut_rhs); +#else + glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP, + mir->cut_rhs); +#endif +#else + { int i; + i = glp_add_rows(pool, 1); + glp_set_row_bnds(pool, i, GLP_UP, 0, mir->cut_rhs); + glp_set_mat_row(pool, i, len, ind, val); + } +#endif + xfree(ind); + xfree(val); + return; +} + +#if 0 /* 29/II-2016 by Chris */ +static int aggregate_row(glp_prob *mip, glp_mir *mir) +#else +static int aggregate_row(glp_prob *mip, glp_mir *mir, SPV *v) +#endif +{ /* try to aggregate another row */ + int m = mir->m; + int n = mir->n; + GLPAIJ *aij; +#if 0 /* 29/II-2016 by Chris */ + SPV *v; +#endif + int ii, j, jj, k, kk, kappa = 0, ret = 0; + double d1, d2, d, d_max = 0.0; + /* choose appropriate structural variable in the aggregated row + to be substituted */ + for (j = 1; j <= mir->agg_vec->nnz; j++) + { k = mir->agg_vec->ind[j]; + xassert(1 <= k && k <= m+n); + if (k <= m) continue; /* skip auxiliary var */ + if (mir->isint[k]) continue; /* skip integer var */ + if (fabs(mir->agg_vec->val[j]) < 0.001) continue; + /* 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); + /* d = min(d1, d2) */ + d = (d1 <= d2 ? d1 : d2); + xassert(d != DBL_MAX); + /* should not be close to corresponding bound */ + if (d < 0.001) continue; + if (d_max < d) d_max = d, kappa = k; + } + if (kappa == 0) + { /* nothing chosen */ + ret = 1; + goto done; + } + /* x[kappa] has been chosen */ + xassert(m+1 <= kappa && kappa <= m+n); + xassert(!mir->isint[kappa]); + /* find another row, which have not been used yet, to eliminate + x[kappa] from the aggregated row */ +#if 0 /* 29/II-2016 by Chris */ + for (ii = 1; ii <= m; ii++) + { if (mir->skip[ii]) continue; + for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next) + if (aij->col->j == kappa - m) break; + if (aij != NULL && fabs(aij->val) >= 0.001) break; +#else + ii = 0; + for (aij = mip->col[kappa - m]->ptr; aij != NULL; + aij = aij->c_next) + { if (aij->row->i > m) continue; + if (mir->skip[aij->row->i]) continue; + if (fabs(aij->val) >= 0.001) + { ii = aij->row->i; + break; + } +#endif + } +#if 0 /* 29/II-2016 by Chris */ + if (ii > m) +#else + if (ii == 0) +#endif + { /* nothing found */ + ret = 2; + goto done; + } + /* row ii has been found; include it in the aggregated list */ + mir->agg_cnt++; + xassert(mir->agg_cnt <= MAXAGGR); + mir->agg_row[mir->agg_cnt] = ii; + mir->skip[ii] = 2; + /* v := new row */ +#if 0 /* 29/II-2016 by Chris */ + v = ios_create_vec(m+n); +#else + spv_clear_vec(v); +#endif + spv_set_vj(v, ii, 1.0); + for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next) + spv_set_vj(v, m + aij->col->j, - aij->val); +#if MIR_DEBUG + spv_check_vec(v); +#endif + /* perform gaussian elimination to remove x[kappa] */ + j = mir->agg_vec->pos[kappa]; + xassert(j != 0); + jj = v->pos[kappa]; + xassert(jj != 0); + spv_linear_comb(mir->agg_vec, + - mir->agg_vec->val[j] / v->val[jj], v); +#if 0 /* 29/II-2016 by Chris */ + ios_delete_vec(v); +#endif + spv_set_vj(mir->agg_vec, kappa, 0.0); +#if MIR_DEBUG + spv_check_vec(mir->agg_vec); +#endif +done: return ret; +} + +int glp_mir_gen(glp_prob *mip, glp_mir *mir, glp_prob *pool) +{ /* main routine to generate MIR cuts */ + int m = mir->m; + int n = mir->n; + int i, nnn = 0; + double r_best; +#if 1 /* 29/II-2016 by Chris */ + SPV *work; +#endif + xassert(mip->m >= m); + xassert(mip->n == n); + /* obtain current point */ + get_current_point(mip, mir); +#if MIR_DEBUG + /* check current point */ + check_current_point(mir); +#endif + /* reset bound substitution flags */ + memset(&mir->subst[1], '?', m+n); +#if 1 /* 29/II-2016 by Chris */ + work = spv_create_vec(m+n); +#endif + /* try to generate a set of violated MIR cuts */ + for (i = 1; i <= m; i++) + { if (mir->skip[i]) continue; + /* use original i-th row as initial aggregated constraint */ + initial_agg_row(mip, mir, i); +loop: ; +#if MIR_DEBUG + /* check aggregated row */ + check_agg_row(mir); +#endif + /* substitute fixed variables into aggregated constraint */ + subst_fixed_vars(mir); +#if MIR_DEBUG + /* check aggregated row */ + check_agg_row(mir); +#endif +#if MIR_DEBUG + /* check bound substitution flags */ + { int k; + for (k = 1; k <= m+n; k++) + xassert(mir->subst[k] == '?'); + } +#endif + /* apply bound substitution heuristic */ + bound_subst_heur(mir); + /* substitute bounds and build modified constraint */ + build_mod_row(mir); +#if MIR_DEBUG + /* check modified row */ + check_mod_row(mir); +#endif + /* try to generate violated c-MIR cut for modified row */ + r_best = generate(mir); + if (r_best > 0.0) + { /* success */ +#if MIR_DEBUG + /* check raw cut before back bound substitution */ + check_raw_cut(mir, r_best); +#endif + /* back substitution of original bounds */ + back_subst(mir); +#if MIR_DEBUG + /* check the cut after back bound substitution */ + check_cut_row(mir, r_best); +#endif + /* final substitution to eliminate auxiliary variables */ + subst_aux_vars(mip, mir); +#if MIR_DEBUG + /* check the cut after elimination of auxiliaries */ + check_cut_row(mir, r_best); +#endif + /* add constructed cut inequality to the cut pool */ + add_cut(mir, pool), nnn++; + } + /* reset bound substitution flags */ + { int j, k; + for (j = 1; j <= mir->mod_vec->nnz; j++) + { k = mir->mod_vec->ind[j]; + xassert(1 <= k && k <= m+n); + xassert(mir->subst[k] != '?'); + mir->subst[k] = '?'; + } + } + if (r_best == 0.0) + { /* failure */ + if (mir->agg_cnt < MAXAGGR) + { /* try to aggregate another row */ +#if 0 /* 29/II-2016 by Chris */ + if (aggregate_row(mip, mir) == 0) goto loop; +#else + if (aggregate_row(mip, mir, work) == 0) goto loop; +#endif + } + } + /* unmark rows used in the aggregated constraint */ + { int k, ii; + for (k = 1; k <= mir->agg_cnt; k++) + { ii = mir->agg_row[k]; + xassert(1 <= ii && ii <= m); + xassert(mir->skip[ii] == 2); + mir->skip[ii] = 0; + } + } + } +#if 1 /* 29/II-2016 by Chris */ + spv_delete_vec(work); +#endif + return nnn; +} + +/*********************************************************************** +* NAME +* +* glp_mir_free - delete MIR cut generator workspace +* +* SYNOPSIS +* +* void glp_mir_free(glp_mir *mir); +* +* DESCRIPTION +* +* This routine deletes the MIR cut generator workspace and frees all +* the memory allocated to it. */ + +void glp_mir_free(glp_mir *mir) +{ xfree(mir->skip); + xfree(mir->isint); + xfree(mir->lb); + xfree(mir->vlb); + xfree(mir->ub); + xfree(mir->vub); + xfree(mir->x); + xfree(mir->agg_row); + spv_delete_vec(mir->agg_vec); + xfree(mir->subst); + spv_delete_vec(mir->mod_vec); + spv_delete_vec(mir->cut_vec); + xfree(mir); + return; +} + +/* eof */ diff --git a/test/monniaux/glpk-4.65/src/intopt/spv.c b/test/monniaux/glpk-4.65/src/intopt/spv.c new file mode 100644 index 00000000..502f3cd0 --- /dev/null +++ b/test/monniaux/glpk-4.65/src/intopt/spv.c @@ -0,0 +1,303 @@ +/* spv.c (operations on sparse vectors) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2007-2017 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 "env.h" +#include "spv.h" + +/*********************************************************************** +* NAME +* +* spv_create_vec - create sparse vector +* +* SYNOPSIS +* +* #include "glpios.h" +* SPV *spv_create_vec(int n); +* +* DESCRIPTION +* +* The routine spv_create_vec creates a sparse vector of dimension n, +* which initially is a null vector. +* +* RETURNS +* +* The routine returns a pointer to the vector created. */ + +SPV *spv_create_vec(int n) +{ SPV *v; + xassert(n >= 0); + v = xmalloc(sizeof(SPV)); + v->n = n; + v->nnz = 0; + v->pos = xcalloc(1+n, sizeof(int)); + memset(&v->pos[1], 0, n * sizeof(int)); + v->ind = xcalloc(1+n, sizeof(int)); + v->val = xcalloc(1+n, sizeof(double)); + return v; +} + +/*********************************************************************** +* NAME +* +* spv_check_vec - check that sparse vector has correct representation +* +* SYNOPSIS +* +* #include "glpios.h" +* void spv_check_vec(SPV *v); +* +* DESCRIPTION +* +* The routine spv_check_vec checks that a sparse vector specified by +* the parameter v has correct representation. +* +* NOTE +* +* Complexity of this operation is O(n). */ + +void spv_check_vec(SPV *v) +{ int j, k, nnz; + xassert(v->n >= 0); + nnz = 0; + for (j = v->n; j >= 1; j--) + { k = v->pos[j]; + xassert(0 <= k && k <= v->nnz); + if (k != 0) + { xassert(v->ind[k] == j); + nnz++; + } + } + xassert(v->nnz == nnz); + return; +} + +/*********************************************************************** +* NAME +* +* spv_get_vj - retrieve component of sparse vector +* +* SYNOPSIS +* +* #include "glpios.h" +* double spv_get_vj(SPV *v, int j); +* +* RETURNS +* +* The routine spv_get_vj returns j-th component of a sparse vector +* specified by the parameter v. */ + +double spv_get_vj(SPV *v, int j) +{ int k; + xassert(1 <= j && j <= v->n); + k = v->pos[j]; + xassert(0 <= k && k <= v->nnz); + return (k == 0 ? 0.0 : v->val[k]); +} + +/*********************************************************************** +* NAME +* +* spv_set_vj - set/change component of sparse vector +* +* SYNOPSIS +* +* #include "glpios.h" +* void spv_set_vj(SPV *v, int j, double val); +* +* DESCRIPTION +* +* The routine spv_set_vj assigns val to j-th component of a sparse +* vector specified by the parameter v. */ + +void spv_set_vj(SPV *v, int j, double val) +{ int k; + xassert(1 <= j && j <= v->n); + k = v->pos[j]; + if (val == 0.0) + { if (k != 0) + { /* remove j-th component */ + v->pos[j] = 0; + if (k < v->nnz) + { v->pos[v->ind[v->nnz]] = k; + v->ind[k] = v->ind[v->nnz]; + v->val[k] = v->val[v->nnz]; + } + v->nnz--; + } + } + else + { if (k == 0) + { /* create j-th component */ + k = ++(v->nnz); + v->pos[j] = k; + v->ind[k] = j; + } + v->val[k] = val; + } + return; +} + +/*********************************************************************** +* NAME +* +* spv_clear_vec - set all components of sparse vector to zero +* +* SYNOPSIS +* +* #include "glpios.h" +* void spv_clear_vec(SPV *v); +* +* DESCRIPTION +* +* The routine spv_clear_vec sets all components of a sparse vector +* specified by the parameter v to zero. */ + +void spv_clear_vec(SPV *v) +{ int k; + for (k = 1; k <= v->nnz; k++) + v->pos[v->ind[k]] = 0; + v->nnz = 0; + return; +} + +/*********************************************************************** +* NAME +* +* spv_clean_vec - remove zero or small components from sparse vector +* +* SYNOPSIS +* +* #include "glpios.h" +* void spv_clean_vec(SPV *v, double eps); +* +* DESCRIPTION +* +* The routine spv_clean_vec removes zero components and components +* whose magnitude is less than eps from a sparse vector specified by +* the parameter v. If eps is 0.0, only zero components are removed. */ + +void spv_clean_vec(SPV *v, double eps) +{ int k, nnz; + nnz = 0; + for (k = 1; k <= v->nnz; k++) + { if (fabs(v->val[k]) == 0.0 || fabs(v->val[k]) < eps) + { /* remove component */ + v->pos[v->ind[k]] = 0; + } + else + { /* keep component */ + nnz++; + v->pos[v->ind[k]] = nnz; + v->ind[nnz] = v->ind[k]; + v->val[nnz] = v->val[k]; + } + } + v->nnz = nnz; + return; +} + +/*********************************************************************** +* NAME +* +* spv_copy_vec - copy sparse vector (x := y) +* +* SYNOPSIS +* +* #include "glpios.h" +* void spv_copy_vec(SPV *x, SPV *y); +* +* DESCRIPTION +* +* The routine spv_copy_vec copies a sparse vector specified by the +* parameter y to a sparse vector specified by the parameter x. */ + +void spv_copy_vec(SPV *x, SPV *y) +{ int j; + xassert(x != y); + xassert(x->n == y->n); + spv_clear_vec(x); + x->nnz = y->nnz; + memcpy(&x->ind[1], &y->ind[1], x->nnz * sizeof(int)); + memcpy(&x->val[1], &y->val[1], x->nnz * sizeof(double)); + for (j = 1; j <= x->nnz; j++) + x->pos[x->ind[j]] = j; + return; +} + +/*********************************************************************** +* NAME +* +* spv_linear_comb - compute linear combination (x := x + a * y) +* +* SYNOPSIS +* +* #include "glpios.h" +* void spv_linear_comb(SPV *x, double a, SPV *y); +* +* DESCRIPTION +* +* The routine spv_linear_comb computes the linear combination +* +* x := x + a * y, +* +* where x and y are sparse vectors, a is a scalar. */ + +void spv_linear_comb(SPV *x, double a, SPV *y) +{ int j, k; + double xj, yj; + xassert(x != y); + xassert(x->n == y->n); + for (k = 1; k <= y->nnz; k++) + { j = y->ind[k]; + xj = spv_get_vj(x, j); + yj = y->val[k]; + spv_set_vj(x, j, xj + a * yj); + } + return; +} + +/*********************************************************************** +* NAME +* +* spv_delete_vec - delete sparse vector +* +* SYNOPSIS +* +* #include "glpios.h" +* void spv_delete_vec(SPV *v); +* +* DESCRIPTION +* +* The routine spv_delete_vec deletes a sparse vector specified by the +* parameter v freeing all the memory allocated to this object. */ + +void spv_delete_vec(SPV *v) +{ /* delete sparse vector */ + xfree(v->pos); + xfree(v->ind); + xfree(v->val); + xfree(v); + return; +} + +/* eof */ diff --git a/test/monniaux/glpk-4.65/src/intopt/spv.h b/test/monniaux/glpk-4.65/src/intopt/spv.h new file mode 100644 index 00000000..d7d4699f --- /dev/null +++ b/test/monniaux/glpk-4.65/src/intopt/spv.h @@ -0,0 +1,83 @@ +/* spv.h (operations on sparse vectors) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2007-2017 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/>. +***********************************************************************/ + +#ifndef SPV_H +#define SPV_H + +typedef struct SPV SPV; + +struct SPV +{ /* sparse vector v = (v[j]) */ + int n; + /* dimension, n >= 0 */ + int nnz; + /* number of non-zero components, 0 <= nnz <= n */ + int *pos; /* int pos[1+n]; */ + /* pos[j] = k, 1 <= j <= n, is position of (non-zero) v[j] in the + * arrays ind and val, where 1 <= k <= nnz; pos[j] = 0 means that + * v[j] is structural zero */ + int *ind; /* int ind[1+n]; */ + /* ind[k] = j, 1 <= k <= nnz, is index of v[j] */ + double *val; /* double val[1+n]; */ + /* val[k], 1 <= k <= nnz, is a numeric value of v[j] */ +}; + +#define spv_create_vec _glp_spv_create_vec +SPV *spv_create_vec(int n); +/* create sparse vector */ + +#define spv_check_vec _glp_spv_check_vec +void spv_check_vec(SPV *v); +/* check that sparse vector has correct representation */ + +#define spv_get_vj _glp_spv_get_vj +double spv_get_vj(SPV *v, int j); +/* retrieve component of sparse vector */ + +#define spv_set_vj _glp_spv_set_vj +void spv_set_vj(SPV *v, int j, double val); +/* set/change component of sparse vector */ + +#define spv_clear_vec _glp_spv_clear_vec +void spv_clear_vec(SPV *v); +/* set all components of sparse vector to zero */ + +#define spv_clean_vec _glp_spv_clean_vec +void spv_clean_vec(SPV *v, double eps); +/* remove zero or small components from sparse vector */ + +#define spv_copy_vec _glp_spv_copy_vec +void spv_copy_vec(SPV *x, SPV *y); +/* copy sparse vector (x := y) */ + +#define spv_linear_comb _glp_spv_linear_comb +void spv_linear_comb(SPV *x, double a, SPV *y); +/* compute linear combination (x := x + a * y) */ + +#define spv_delete_vec _glp_spv_delete_vec +void spv_delete_vec(SPV *v); +/* delete sparse vector */ + +#endif + +/* eof */ |