alpar@9: /* glpnpp04.c */ alpar@9: alpar@9: /*********************************************************************** alpar@9: * This code is part of GLPK (GNU Linear Programming Kit). alpar@9: * alpar@9: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@9: * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, alpar@9: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@9: * E-mail: . alpar@9: * alpar@9: * GLPK is free software: you can redistribute it and/or modify it alpar@9: * under the terms of the GNU General Public License as published by alpar@9: * the Free Software Foundation, either version 3 of the License, or alpar@9: * (at your option) any later version. alpar@9: * alpar@9: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@9: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@9: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@9: * License for more details. alpar@9: * alpar@9: * You should have received a copy of the GNU General Public License alpar@9: * along with GLPK. If not, see . alpar@9: ***********************************************************************/ alpar@9: alpar@9: #include "glpnpp.h" alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_binarize_prob - binarize MIP problem alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_binarize_prob(NPP *npp); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_binarize_prob replaces in the original MIP problem alpar@9: * every integer variable: alpar@9: * alpar@9: * l[q] <= x[q] <= u[q], (1) alpar@9: * alpar@9: * where l[q] < u[q], by an equivalent sum of binary variables. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * The routine returns the number of integer variables for which the alpar@9: * transformation failed, because u[q] - l[q] > d_max. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * If variable x[q] has non-zero lower bound, it is first processed alpar@9: * with the routine npp_lbnd_col. Thus, we can assume that: alpar@9: * alpar@9: * 0 <= x[q] <= u[q]. (2) alpar@9: * alpar@9: * If u[q] = 1, variable x[q] is already binary, so further processing alpar@9: * is not needed. Let, therefore, that 2 <= u[q] <= d_max, and n be a alpar@9: * smallest integer such that u[q] <= 2^n - 1 (n >= 2, since u[q] >= 2). alpar@9: * Then variable x[q] can be replaced by the following sum: alpar@9: * alpar@9: * n-1 alpar@9: * x[q] = sum 2^k x[k], (3) alpar@9: * k=0 alpar@9: * alpar@9: * where x[k] are binary columns (variables). If u[q] < 2^n - 1, the alpar@9: * following additional inequality constraint must be also included in alpar@9: * the transformed problem: alpar@9: * alpar@9: * n-1 alpar@9: * sum 2^k x[k] <= u[q]. (4) alpar@9: * k=0 alpar@9: * alpar@9: * Note: Assuming that in the transformed problem x[q] becomes binary alpar@9: * variable x[0], this transformation causes new n-1 binary variables alpar@9: * to appear. alpar@9: * alpar@9: * Substituting x[q] from (3) to the objective row gives: alpar@9: * alpar@9: * z = sum c[j] x[j] + c[0] = alpar@9: * j alpar@9: * alpar@9: * = sum c[j] x[j] + c[q] x[q] + c[0] = alpar@9: * j!=q alpar@9: * n-1 alpar@9: * = sum c[j] x[j] + c[q] sum 2^k x[k] + c[0] = alpar@9: * j!=q k=0 alpar@9: * n-1 alpar@9: * = sum c[j] x[j] + sum c[k] x[k] + c[0], alpar@9: * j!=q k=0 alpar@9: * alpar@9: * where: alpar@9: * alpar@9: * c[k] = 2^k c[q], k = 0, ..., n-1. (5) alpar@9: * alpar@9: * And substituting x[q] from (3) to i-th constraint row i gives: alpar@9: * alpar@9: * L[i] <= sum a[i,j] x[j] <= U[i] ==> alpar@9: * j alpar@9: * alpar@9: * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==> alpar@9: * j!=q alpar@9: * n-1 alpar@9: * L[i] <= sum a[i,j] x[j] + a[i,q] sum 2^k x[k] <= U[i] ==> alpar@9: * j!=q k=0 alpar@9: * n-1 alpar@9: * L[i] <= sum a[i,j] x[j] + sum a[i,k] x[k] <= U[i], alpar@9: * j!=q k=0 alpar@9: * alpar@9: * where: alpar@9: * alpar@9: * a[i,k] = 2^k a[i,q], k = 0, ..., n-1. (6) alpar@9: * alpar@9: * RECOVERING SOLUTION alpar@9: * alpar@9: * Value of variable x[q] is computed with formula (3). */ alpar@9: alpar@9: struct binarize alpar@9: { int q; alpar@9: /* column reference number for x[q] = x[0] */ alpar@9: int j; alpar@9: /* column reference number for x[1]; x[2] has reference number alpar@9: j+1, x[3] - j+2, etc. */ alpar@9: int n; alpar@9: /* total number of binary variables, n >= 2 */ alpar@9: }; alpar@9: alpar@9: static int rcv_binarize_prob(NPP *npp, void *info); alpar@9: alpar@9: int npp_binarize_prob(NPP *npp) alpar@9: { /* binarize MIP problem */ alpar@9: struct binarize *info; alpar@9: NPPROW *row; alpar@9: NPPCOL *col, *bin; alpar@9: NPPAIJ *aij; alpar@9: int u, n, k, temp, nfails, nvars, nbins, nrows; alpar@9: /* new variables will be added to the end of the column list, so alpar@9: we go from the end to beginning of the column list */ alpar@9: nfails = nvars = nbins = nrows = 0; alpar@9: for (col = npp->c_tail; col != NULL; col = col->prev) alpar@9: { /* skip continuous variable */ alpar@9: if (!col->is_int) continue; alpar@9: /* skip fixed variable */ alpar@9: if (col->lb == col->ub) continue; alpar@9: /* skip binary variable */ alpar@9: if (col->lb == 0.0 && col->ub == 1.0) continue; alpar@9: /* check if the transformation is applicable */ alpar@9: if (col->lb < -1e6 || col->ub > +1e6 || alpar@9: col->ub - col->lb > 4095.0) alpar@9: { /* unfortunately, not */ alpar@9: nfails++; alpar@9: continue; alpar@9: } alpar@9: /* process integer non-binary variable x[q] */ alpar@9: nvars++; alpar@9: /* make x[q] non-negative, if its lower bound is non-zero */ alpar@9: if (col->lb != 0.0) alpar@9: npp_lbnd_col(npp, col); alpar@9: /* now 0 <= x[q] <= u[q] */ alpar@9: xassert(col->lb == 0.0); alpar@9: u = (int)col->ub; alpar@9: xassert(col->ub == (double)u); alpar@9: /* if x[q] is binary, further processing is not needed */ alpar@9: if (u == 1) continue; alpar@9: /* determine smallest n such that u <= 2^n - 1 (thus, n is the alpar@9: number of binary variables needed) */ alpar@9: n = 2, temp = 4; alpar@9: while (u >= temp) alpar@9: n++, temp += temp; alpar@9: nbins += n; alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_binarize_prob, sizeof(struct binarize)); alpar@9: info->q = col->j; alpar@9: info->j = 0; /* will be set below */ alpar@9: info->n = n; alpar@9: /* if u < 2^n - 1, we need one additional row for (4) */ alpar@9: if (u < temp - 1) alpar@9: { row = npp_add_row(npp), nrows++; alpar@9: row->lb = -DBL_MAX, row->ub = u; alpar@9: } alpar@9: else alpar@9: row = NULL; alpar@9: /* in the transformed problem variable x[q] becomes binary alpar@9: variable x[0], so its objective and constraint coefficients alpar@9: are not changed */ alpar@9: col->ub = 1.0; alpar@9: /* include x[0] into constraint (4) */ alpar@9: if (row != NULL) alpar@9: npp_add_aij(npp, row, col, 1.0); alpar@9: /* add other binary variables x[1], ..., x[n-1] */ alpar@9: for (k = 1, temp = 2; k < n; k++, temp += temp) alpar@9: { /* add new binary variable x[k] */ alpar@9: bin = npp_add_col(npp); alpar@9: bin->is_int = 1; alpar@9: bin->lb = 0.0, bin->ub = 1.0; alpar@9: bin->coef = (double)temp * col->coef; alpar@9: /* store column reference number for x[1] */ alpar@9: if (info->j == 0) alpar@9: info->j = bin->j; alpar@9: else alpar@9: xassert(info->j + (k-1) == bin->j); alpar@9: /* duplicate constraint coefficients for x[k]; this also alpar@9: automatically includes x[k] into constraint (4) */ alpar@9: for (aij = col->ptr; aij != NULL; aij = aij->c_next) alpar@9: npp_add_aij(npp, aij->row, bin, (double)temp * aij->val); alpar@9: } alpar@9: } alpar@9: if (nvars > 0) alpar@9: xprintf("%d integer variable(s) were replaced by %d binary one" alpar@9: "s\n", nvars, nbins); alpar@9: if (nrows > 0) alpar@9: xprintf("%d row(s) were added due to binarization\n", nrows); alpar@9: if (nfails > 0) alpar@9: xprintf("Binarization failed for %d integer variable(s)\n", alpar@9: nfails); alpar@9: return nfails; alpar@9: } alpar@9: alpar@9: static int rcv_binarize_prob(NPP *npp, void *_info) alpar@9: { /* recovery binarized variable */ alpar@9: struct binarize *info = _info; alpar@9: int k, temp; alpar@9: double sum; alpar@9: /* compute value of x[q]; see formula (3) */ alpar@9: sum = npp->c_value[info->q]; alpar@9: for (k = 1, temp = 2; k < info->n; k++, temp += temp) alpar@9: sum += (double)temp * npp->c_value[info->j + (k-1)]; alpar@9: npp->c_value[info->q] = sum; alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /**********************************************************************/ alpar@9: alpar@9: struct elem alpar@9: { /* linear form element a[j] x[j] */ alpar@9: double aj; alpar@9: /* non-zero coefficient value */ alpar@9: NPPCOL *xj; alpar@9: /* pointer to variable (column) */ alpar@9: struct elem *next; alpar@9: /* pointer to another term */ alpar@9: }; alpar@9: alpar@9: static struct elem *copy_form(NPP *npp, NPPROW *row, double s) alpar@9: { /* copy linear form */ alpar@9: NPPAIJ *aij; alpar@9: struct elem *ptr, *e; alpar@9: ptr = NULL; alpar@9: for (aij = row->ptr; aij != NULL; aij = aij->r_next) alpar@9: { e = dmp_get_atom(npp->pool, sizeof(struct elem)); alpar@9: e->aj = s * aij->val; alpar@9: e->xj = aij->col; alpar@9: e->next = ptr; alpar@9: ptr = e; alpar@9: } alpar@9: return ptr; alpar@9: } alpar@9: alpar@9: static void drop_form(NPP *npp, struct elem *ptr) alpar@9: { /* drop linear form */ alpar@9: struct elem *e; alpar@9: while (ptr != NULL) alpar@9: { e = ptr; alpar@9: ptr = e->next; alpar@9: dmp_free_atom(npp->pool, e, sizeof(struct elem)); alpar@9: } alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_is_packing - test if constraint is packing inequality alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_is_packing(NPP *npp, NPPROW *row); alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * If the specified row (constraint) is packing inequality (see below), alpar@9: * the routine npp_is_packing returns non-zero. Otherwise, it returns alpar@9: * zero. alpar@9: * alpar@9: * PACKING INEQUALITIES alpar@9: * alpar@9: * In canonical format the packing inequality is the following: alpar@9: * alpar@9: * sum x[j] <= 1, (1) alpar@9: * j in J alpar@9: * alpar@9: * where all variables x[j] are binary. This inequality expresses the alpar@9: * condition that in any integer feasible solution at most one variable alpar@9: * from set J can take non-zero (unity) value while other variables alpar@9: * must be equal to zero. W.l.o.g. it is assumed that |J| >= 2, because alpar@9: * if J is empty or |J| = 1, the inequality (1) is redundant. alpar@9: * alpar@9: * In general case the packing inequality may include original variables alpar@9: * x[j] as well as their complements x~[j]: alpar@9: * alpar@9: * sum x[j] + sum x~[j] <= 1, (2) alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * where Jp and Jn are not intersected. Therefore, using substitution alpar@9: * x~[j] = 1 - x[j] gives the packing inequality in generalized format: alpar@9: * alpar@9: * sum x[j] - sum x[j] <= 1 - |Jn|. (3) alpar@9: * j in Jp j in Jn */ alpar@9: alpar@9: int npp_is_packing(NPP *npp, NPPROW *row) alpar@9: { /* test if constraint is packing inequality */ alpar@9: NPPCOL *col; alpar@9: NPPAIJ *aij; alpar@9: int b; alpar@9: xassert(npp == npp); alpar@9: if (!(row->lb == -DBL_MAX && row->ub != +DBL_MAX)) alpar@9: return 0; alpar@9: b = 1; alpar@9: for (aij = row->ptr; aij != NULL; aij = aij->r_next) alpar@9: { col = aij->col; alpar@9: if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0)) alpar@9: return 0; alpar@9: if (aij->val == +1.0) alpar@9: ; alpar@9: else if (aij->val == -1.0) alpar@9: b--; alpar@9: else alpar@9: return 0; alpar@9: } alpar@9: if (row->ub != (double)b) return 0; alpar@9: return 1; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_hidden_packing - identify hidden packing inequality alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_hidden_packing(NPP *npp, NPPROW *row); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_hidden_packing processes specified inequality alpar@9: * constraint, which includes only binary variables, and the number of alpar@9: * the variables is not less than two. If the original inequality is alpar@9: * equivalent to a packing inequality, the routine replaces it by this alpar@9: * equivalent inequality. If the original constraint is double-sided alpar@9: * inequality, it is replaced by a pair of single-sided inequalities, alpar@9: * if necessary. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * If the original inequality constraint was replaced by equivalent alpar@9: * packing inequality, the routine npp_hidden_packing returns non-zero. alpar@9: * Otherwise, it returns zero. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * Consider an inequality constraint: alpar@9: * alpar@9: * sum a[j] x[j] <= b, (1) alpar@9: * j in J alpar@9: * alpar@9: * where all variables x[j] are binary, and |J| >= 2. (In case of '>=' alpar@9: * inequality it can be transformed to '<=' format by multiplying both alpar@9: * its sides by -1.) alpar@9: * alpar@9: * Let Jp = {j: a[j] > 0}, Jn = {j: a[j] < 0}. Performing substitution alpar@9: * x[j] = 1 - x~[j] for all j in Jn, we have: alpar@9: * alpar@9: * sum a[j] x[j] <= b ==> alpar@9: * j in J alpar@9: * alpar@9: * sum a[j] x[j] + sum a[j] x[j] <= b ==> alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * sum a[j] x[j] + sum a[j] (1 - x~[j]) <= b ==> alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * sum a[j] x[j] - sum a[j] x~[j] <= b - sum a[j]. alpar@9: * j in Jp j in Jn j in Jn alpar@9: * alpar@9: * Thus, meaning the transformation above, we can assume that in alpar@9: * inequality (1) all coefficients a[j] are positive. Moreover, we can alpar@9: * assume that a[j] <= b. In fact, let a[j] > b; then the following alpar@9: * three cases are possible: alpar@9: * alpar@9: * 1) b < 0. In this case inequality (1) is infeasible, so the problem alpar@9: * has no feasible solution (see the routine npp_analyze_row); alpar@9: * alpar@9: * 2) b = 0. In this case inequality (1) is a forcing inequality on its alpar@9: * upper bound (see the routine npp_forcing row), from which it alpar@9: * follows that all variables x[j] should be fixed at zero; alpar@9: * alpar@9: * 3) b > 0. In this case inequality (1) defines an implied zero upper alpar@9: * bound for variable x[j] (see the routine npp_implied_bounds), from alpar@9: * which it follows that x[j] should be fixed at zero. alpar@9: * alpar@9: * It is assumed that all three cases listed above have been recognized alpar@9: * by the routine npp_process_prob, which performs basic MIP processing alpar@9: * prior to a call the routine npp_hidden_packing. So, if one of these alpar@9: * cases occurs, we should just skip processing such constraint. alpar@9: * alpar@9: * Thus, let 0 < a[j] <= b. Then it is obvious that constraint (1) is alpar@9: * equivalent to packing inquality only if: alpar@9: * alpar@9: * a[j] + a[k] > b + eps (2) alpar@9: * alpar@9: * for all j, k in J, j != k, where eps is an absolute tolerance for alpar@9: * row (linear form) value. Checking the condition (2) for all j and k, alpar@9: * j != k, requires time O(|J|^2). However, this time can be reduced to alpar@9: * O(|J|), if use minimal a[j] and a[k], in which case it is sufficient alpar@9: * to check the condition (2) only once. alpar@9: * alpar@9: * Once the original inequality (1) is replaced by equivalent packing alpar@9: * inequality, we need to perform back substitution x~[j] = 1 - x[j] for alpar@9: * all j in Jn (see above). alpar@9: * alpar@9: * RECOVERING SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: static int hidden_packing(NPP *npp, struct elem *ptr, double *_b) alpar@9: { /* process inequality constraint: sum a[j] x[j] <= b; alpar@9: 0 - specified row is NOT hidden packing inequality; alpar@9: 1 - specified row is packing inequality; alpar@9: 2 - specified row is hidden packing inequality. */ alpar@9: struct elem *e, *ej, *ek; alpar@9: int neg; alpar@9: double b = *_b, eps; alpar@9: xassert(npp == npp); alpar@9: /* a[j] must be non-zero, x[j] must be binary, for all j in J */ alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: { xassert(e->aj != 0.0); alpar@9: xassert(e->xj->is_int); alpar@9: xassert(e->xj->lb == 0.0 && e->xj->ub == 1.0); alpar@9: } alpar@9: /* check if the specified inequality constraint already has the alpar@9: form of packing inequality */ alpar@9: neg = 0; /* neg is |Jn| */ alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: { if (e->aj == +1.0) alpar@9: ; alpar@9: else if (e->aj == -1.0) alpar@9: neg++; alpar@9: else alpar@9: break; alpar@9: } alpar@9: if (e == NULL) alpar@9: { /* all coefficients a[j] are +1 or -1; check rhs b */ alpar@9: if (b == (double)(1 - neg)) alpar@9: { /* it is packing inequality; no processing is needed */ alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: /* substitute x[j] = 1 - x~[j] for all j in Jn to make all a[j] alpar@9: positive; the result is a~[j] = |a[j]| and new rhs b */ alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: if (e->aj < 0) b -= e->aj; alpar@9: /* now a[j] > 0 for all j in J (actually |a[j]| are used) */ alpar@9: /* if a[j] > b, skip processing--this case must not appear */ alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: if (fabs(e->aj) > b) return 0; alpar@9: /* now 0 < a[j] <= b for all j in J */ alpar@9: /* find two minimal coefficients a[j] and a[k], j != k */ alpar@9: ej = NULL; alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: if (ej == NULL || fabs(ej->aj) > fabs(e->aj)) ej = e; alpar@9: xassert(ej != NULL); alpar@9: ek = NULL; alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: if (e != ej) alpar@9: if (ek == NULL || fabs(ek->aj) > fabs(e->aj)) ek = e; alpar@9: xassert(ek != NULL); alpar@9: /* the specified constraint is equivalent to packing inequality alpar@9: iff a[j] + a[k] > b + eps */ alpar@9: eps = 1e-3 + 1e-6 * fabs(b); alpar@9: if (fabs(ej->aj) + fabs(ek->aj) <= b + eps) return 0; alpar@9: /* perform back substitution x~[j] = 1 - x[j] and construct the alpar@9: final equivalent packing inequality in generalized format */ alpar@9: b = 1.0; alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: { if (e->aj > 0.0) alpar@9: e->aj = +1.0; alpar@9: else /* e->aj < 0.0 */ alpar@9: e->aj = -1.0, b -= 1.0; alpar@9: } alpar@9: *_b = b; alpar@9: return 2; alpar@9: } alpar@9: alpar@9: int npp_hidden_packing(NPP *npp, NPPROW *row) alpar@9: { /* identify hidden packing inequality */ alpar@9: NPPROW *copy; alpar@9: NPPAIJ *aij; alpar@9: struct elem *ptr, *e; alpar@9: int kase, ret, count = 0; alpar@9: double b; alpar@9: /* the row must be inequality constraint */ alpar@9: xassert(row->lb < row->ub); alpar@9: for (kase = 0; kase <= 1; kase++) alpar@9: { if (kase == 0) alpar@9: { /* process row upper bound */ alpar@9: if (row->ub == +DBL_MAX) continue; alpar@9: ptr = copy_form(npp, row, +1.0); alpar@9: b = + row->ub; alpar@9: } alpar@9: else alpar@9: { /* process row lower bound */ alpar@9: if (row->lb == -DBL_MAX) continue; alpar@9: ptr = copy_form(npp, row, -1.0); alpar@9: b = - row->lb; alpar@9: } alpar@9: /* now the inequality has the form "sum a[j] x[j] <= b" */ alpar@9: ret = hidden_packing(npp, ptr, &b); alpar@9: xassert(0 <= ret && ret <= 2); alpar@9: if (kase == 1 && ret == 1 || ret == 2) alpar@9: { /* the original inequality has been identified as hidden alpar@9: packing inequality */ alpar@9: count++; alpar@9: #ifdef GLP_DEBUG alpar@9: xprintf("Original constraint:\n"); alpar@9: for (aij = row->ptr; aij != NULL; aij = aij->r_next) alpar@9: xprintf(" %+g x%d", aij->val, aij->col->j); alpar@9: if (row->lb != -DBL_MAX) xprintf(", >= %g", row->lb); alpar@9: if (row->ub != +DBL_MAX) xprintf(", <= %g", row->ub); alpar@9: xprintf("\n"); alpar@9: xprintf("Equivalent packing inequality:\n"); alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: xprintf(" %sx%d", e->aj > 0.0 ? "+" : "-", e->xj->j); alpar@9: xprintf(", <= %g\n", b); alpar@9: #endif alpar@9: if (row->lb == -DBL_MAX || row->ub == +DBL_MAX) alpar@9: { /* the original row is single-sided inequality; no copy alpar@9: is needed */ alpar@9: copy = NULL; alpar@9: } alpar@9: else alpar@9: { /* the original row is double-sided inequality; we need alpar@9: to create its copy for other bound before replacing it alpar@9: with the equivalent inequality */ alpar@9: copy = npp_add_row(npp); alpar@9: if (kase == 0) alpar@9: { /* the copy is for lower bound */ alpar@9: copy->lb = row->lb, copy->ub = +DBL_MAX; alpar@9: } alpar@9: else alpar@9: { /* the copy is for upper bound */ alpar@9: copy->lb = -DBL_MAX, copy->ub = row->ub; alpar@9: } alpar@9: /* copy original row coefficients */ alpar@9: for (aij = row->ptr; aij != NULL; aij = aij->r_next) alpar@9: npp_add_aij(npp, copy, aij->col, aij->val); alpar@9: } alpar@9: /* replace the original inequality by equivalent one */ alpar@9: npp_erase_row(npp, row); alpar@9: row->lb = -DBL_MAX, row->ub = b; alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: npp_add_aij(npp, row, e->xj, e->aj); alpar@9: /* continue processing lower bound for the copy */ alpar@9: if (copy != NULL) row = copy; alpar@9: } alpar@9: drop_form(npp, ptr); alpar@9: } alpar@9: return count; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_implied_packing - identify implied packing inequality alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_implied_packing(NPP *npp, NPPROW *row, int which, alpar@9: * NPPCOL *var[], char set[]); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_implied_packing processes specified row (constraint) alpar@9: * of general format: alpar@9: * alpar@9: * L <= sum a[j] x[j] <= U. (1) alpar@9: * j alpar@9: * alpar@9: * If which = 0, only lower bound L, which must exist, is considered, alpar@9: * while upper bound U is ignored. Similarly, if which = 1, only upper alpar@9: * bound U, which must exist, is considered, while lower bound L is alpar@9: * ignored. Thus, if the specified row is a double-sided inequality or alpar@9: * equality constraint, this routine should be called twice for both alpar@9: * lower and upper bounds. alpar@9: * alpar@9: * The routine npp_implied_packing attempts to find a non-trivial (i.e. alpar@9: * having not less than two binary variables) packing inequality: alpar@9: * alpar@9: * sum x[j] - sum x[j] <= 1 - |Jn|, (2) alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * which is relaxation of the constraint (1) in the sense that any alpar@9: * solution satisfying to that constraint also satisfies to the packing alpar@9: * inequality (2). If such relaxation exists, the routine stores alpar@9: * pointers to descriptors of corresponding binary variables and their alpar@9: * flags, resp., to locations var[1], var[2], ..., var[len] and set[1], alpar@9: * set[2], ..., set[len], where set[j] = 0 means that j in Jp and alpar@9: * set[j] = 1 means that j in Jn. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * The routine npp_implied_packing returns len, which is the total alpar@9: * number of binary variables in the packing inequality found, len >= 2. alpar@9: * However, if the relaxation does not exist, the routine returns zero. alpar@9: * alpar@9: * ALGORITHM alpar@9: * alpar@9: * If which = 0, the constraint coefficients (1) are multiplied by -1 alpar@9: * and b is assigned -L; if which = 1, the constraint coefficients (1) alpar@9: * are not changed and b is assigned +U. In both cases the specified alpar@9: * constraint gets the following format: alpar@9: * alpar@9: * sum a[j] x[j] <= b. (3) alpar@9: * j alpar@9: * alpar@9: * (Note that (3) is a relaxation of (1), because one of bounds L or U alpar@9: * is ignored.) alpar@9: * alpar@9: * Let J be set of binary variables, Kp be set of non-binary (integer alpar@9: * or continuous) variables with a[j] > 0, and Kn be set of non-binary alpar@9: * variables with a[j] < 0. Then the inequality (3) can be written as alpar@9: * follows: alpar@9: * alpar@9: * sum a[j] x[j] <= b - sum a[j] x[j] - sum a[j] x[j]. (4) alpar@9: * j in J j in Kp j in Kn alpar@9: * alpar@9: * To get rid of non-binary variables we can replace the inequality (4) alpar@9: * by the following relaxed inequality: alpar@9: * alpar@9: * sum a[j] x[j] <= b~, (5) alpar@9: * j in J alpar@9: * alpar@9: * where: alpar@9: * alpar@9: * b~ = sup(b - sum a[j] x[j] - sum a[j] x[j]) = alpar@9: * j in Kp j in Kn alpar@9: * alpar@9: * = b - inf sum a[j] x[j] - inf sum a[j] x[j] = (6) alpar@9: * j in Kp j in Kn alpar@9: * alpar@9: * = b - sum a[j] l[j] - sum a[j] u[j]. alpar@9: * j in Kp j in Kn alpar@9: * alpar@9: * Note that if lower bound l[j] (if j in Kp) or upper bound u[j] alpar@9: * (if j in Kn) of some non-binary variable x[j] does not exist, then alpar@9: * formally b = +oo, in which case further analysis is not performed. alpar@9: * alpar@9: * Let Bp = {j in J: a[j] > 0}, Bn = {j in J: a[j] < 0}. To make all alpar@9: * the inequality coefficients in (5) positive, we replace all x[j] in alpar@9: * Bn by their complementaries, substituting x[j] = 1 - x~[j] for all alpar@9: * j in Bn, that gives: alpar@9: * alpar@9: * sum a[j] x[j] - sum a[j] x~[j] <= b~ - sum a[j]. (7) alpar@9: * j in Bp j in Bn j in Bn alpar@9: * alpar@9: * This inequality is a relaxation of the original constraint (1), and alpar@9: * it is a binary knapsack inequality. Writing it in the standard format alpar@9: * we have: alpar@9: * alpar@9: * sum alfa[j] z[j] <= beta, (8) alpar@9: * j in J alpar@9: * alpar@9: * where: alpar@9: * ( + a[j], if j in Bp, alpar@9: * alfa[j] = < (9) alpar@9: * ( - a[j], if j in Bn, alpar@9: * alpar@9: * ( x[j], if j in Bp, alpar@9: * z[j] = < (10) alpar@9: * ( 1 - x[j], if j in Bn, alpar@9: * alpar@9: * beta = b~ - sum a[j]. (11) alpar@9: * j in Bn alpar@9: * alpar@9: * In the inequality (8) all coefficients are positive, therefore, the alpar@9: * packing relaxation to be found for this inequality is the following: alpar@9: * alpar@9: * sum z[j] <= 1. (12) alpar@9: * j in P alpar@9: * alpar@9: * It is obvious that set P within J, which we would like to find, must alpar@9: * satisfy to the following condition: alpar@9: * alpar@9: * alfa[j] + alfa[k] > beta + eps for all j, k in P, j != k, (13) alpar@9: * alpar@9: * where eps is an absolute tolerance for value of the linear form. alpar@9: * Thus, it is natural to take P = {j: alpha[j] > (beta + eps) / 2}. alpar@9: * Moreover, if in the equality (8) there exist coefficients alfa[k], alpar@9: * for which alfa[k] <= (beta + eps) / 2, but which, nevertheless, alpar@9: * satisfies to the condition (13) for all j in P, *one* corresponding alpar@9: * variable z[k] (having, for example, maximal coefficient alfa[k]) can alpar@9: * be included in set P, that allows increasing the number of binary alpar@9: * variables in (12) by one. alpar@9: * alpar@9: * Once the set P has been built, for the inequality (12) we need to alpar@9: * perform back substitution according to (10) in order to express it alpar@9: * through the original binary variables. As the result of such back alpar@9: * substitution the relaxed packing inequality get its final format (2), alpar@9: * where Jp = J intersect Bp, and Jn = J intersect Bn. */ alpar@9: alpar@9: int npp_implied_packing(NPP *npp, NPPROW *row, int which, alpar@9: NPPCOL *var[], char set[]) alpar@9: { struct elem *ptr, *e, *i, *k; alpar@9: int len = 0; alpar@9: double b, eps; alpar@9: /* build inequality (3) */ alpar@9: if (which == 0) alpar@9: { ptr = copy_form(npp, row, -1.0); alpar@9: xassert(row->lb != -DBL_MAX); alpar@9: b = - row->lb; alpar@9: } alpar@9: else if (which == 1) alpar@9: { ptr = copy_form(npp, row, +1.0); alpar@9: xassert(row->ub != +DBL_MAX); alpar@9: b = + row->ub; alpar@9: } alpar@9: /* remove non-binary variables to build relaxed inequality (5); alpar@9: compute its right-hand side b~ with formula (6) */ alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: { if (!(e->xj->is_int && e->xj->lb == 0.0 && e->xj->ub == 1.0)) alpar@9: { /* x[j] is non-binary variable */ alpar@9: if (e->aj > 0.0) alpar@9: { if (e->xj->lb == -DBL_MAX) goto done; alpar@9: b -= e->aj * e->xj->lb; alpar@9: } alpar@9: else /* e->aj < 0.0 */ alpar@9: { if (e->xj->ub == +DBL_MAX) goto done; alpar@9: b -= e->aj * e->xj->ub; alpar@9: } alpar@9: /* a[j] = 0 means that variable x[j] is removed */ alpar@9: e->aj = 0.0; alpar@9: } alpar@9: } alpar@9: /* substitute x[j] = 1 - x~[j] to build knapsack inequality (8); alpar@9: compute its right-hand side beta with formula (11) */ alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: if (e->aj < 0.0) b -= e->aj; alpar@9: /* if beta is close to zero, the knapsack inequality is either alpar@9: infeasible or forcing inequality; this must never happen, so alpar@9: we skip further analysis */ alpar@9: if (b < 1e-3) goto done; alpar@9: /* build set P as well as sets Jp and Jn, and determine x[k] as alpar@9: explained above in comments to the routine */ alpar@9: eps = 1e-3 + 1e-6 * b; alpar@9: i = k = NULL; alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: { /* note that alfa[j] = |a[j]| */ alpar@9: if (fabs(e->aj) > 0.5 * (b + eps)) alpar@9: { /* alfa[j] > (b + eps) / 2; include x[j] in set P, i.e. in alpar@9: set Jp or Jn */ alpar@9: var[++len] = e->xj; alpar@9: set[len] = (char)(e->aj > 0.0 ? 0 : 1); alpar@9: /* alfa[i] = min alfa[j] over all j included in set P */ alpar@9: if (i == NULL || fabs(i->aj) > fabs(e->aj)) i = e; alpar@9: } alpar@9: else if (fabs(e->aj) >= 1e-3) alpar@9: { /* alfa[k] = max alfa[j] over all j not included in set P; alpar@9: we skip coefficient a[j] if it is close to zero to avoid alpar@9: numerically unreliable results */ alpar@9: if (k == NULL || fabs(k->aj) < fabs(e->aj)) k = e; alpar@9: } alpar@9: } alpar@9: /* if alfa[k] satisfies to condition (13) for all j in P, include alpar@9: x[k] in P */ alpar@9: if (i != NULL && k != NULL && fabs(i->aj) + fabs(k->aj) > b + eps) alpar@9: { var[++len] = k->xj; alpar@9: set[len] = (char)(k->aj > 0.0 ? 0 : 1); alpar@9: } alpar@9: /* trivial packing inequality being redundant must never appear, alpar@9: so we just ignore it */ alpar@9: if (len < 2) len = 0; alpar@9: done: drop_form(npp, ptr); alpar@9: return len; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_is_covering - test if constraint is covering inequality alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_is_covering(NPP *npp, NPPROW *row); alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * If the specified row (constraint) is covering inequality (see below), alpar@9: * the routine npp_is_covering returns non-zero. Otherwise, it returns alpar@9: * zero. alpar@9: * alpar@9: * COVERING INEQUALITIES alpar@9: * alpar@9: * In canonical format the covering inequality is the following: alpar@9: * alpar@9: * sum x[j] >= 1, (1) alpar@9: * j in J alpar@9: * alpar@9: * where all variables x[j] are binary. This inequality expresses the alpar@9: * condition that in any integer feasible solution variables in set J alpar@9: * cannot be all equal to zero at the same time, i.e. at least one alpar@9: * variable must take non-zero (unity) value. W.l.o.g. it is assumed alpar@9: * that |J| >= 2, because if J is empty, the inequality (1) is alpar@9: * infeasible, and if |J| = 1, the inequality (1) is a forcing row. alpar@9: * alpar@9: * In general case the covering inequality may include original alpar@9: * variables x[j] as well as their complements x~[j]: alpar@9: * alpar@9: * sum x[j] + sum x~[j] >= 1, (2) alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * where Jp and Jn are not intersected. Therefore, using substitution alpar@9: * x~[j] = 1 - x[j] gives the packing inequality in generalized format: alpar@9: * alpar@9: * sum x[j] - sum x[j] >= 1 - |Jn|. (3) alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * (May note that the inequality (3) cuts off infeasible solutions, alpar@9: * where x[j] = 0 for all j in Jp and x[j] = 1 for all j in Jn.) alpar@9: * alpar@9: * NOTE: If |J| = 2, the inequality (3) is equivalent to packing alpar@9: * inequality (see the routine npp_is_packing). */ alpar@9: alpar@9: int npp_is_covering(NPP *npp, NPPROW *row) alpar@9: { /* test if constraint is covering inequality */ alpar@9: NPPCOL *col; alpar@9: NPPAIJ *aij; alpar@9: int b; alpar@9: xassert(npp == npp); alpar@9: if (!(row->lb != -DBL_MAX && row->ub == +DBL_MAX)) alpar@9: return 0; alpar@9: b = 1; alpar@9: for (aij = row->ptr; aij != NULL; aij = aij->r_next) alpar@9: { col = aij->col; alpar@9: if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0)) alpar@9: return 0; alpar@9: if (aij->val == +1.0) alpar@9: ; alpar@9: else if (aij->val == -1.0) alpar@9: b--; alpar@9: else alpar@9: return 0; alpar@9: } alpar@9: if (row->lb != (double)b) return 0; alpar@9: return 1; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_hidden_covering - identify hidden covering inequality alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_hidden_covering(NPP *npp, NPPROW *row); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_hidden_covering processes specified inequality alpar@9: * constraint, which includes only binary variables, and the number of alpar@9: * the variables is not less than three. If the original inequality is alpar@9: * equivalent to a covering inequality (see below), the routine alpar@9: * replaces it by the equivalent inequality. If the original constraint alpar@9: * is double-sided inequality, it is replaced by a pair of single-sided alpar@9: * inequalities, if necessary. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * If the original inequality constraint was replaced by equivalent alpar@9: * covering inequality, the routine npp_hidden_covering returns alpar@9: * non-zero. Otherwise, it returns zero. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * Consider an inequality constraint: alpar@9: * alpar@9: * sum a[j] x[j] >= b, (1) alpar@9: * j in J alpar@9: * alpar@9: * where all variables x[j] are binary, and |J| >= 3. (In case of '<=' alpar@9: * inequality it can be transformed to '>=' format by multiplying both alpar@9: * its sides by -1.) alpar@9: * alpar@9: * Let Jp = {j: a[j] > 0}, Jn = {j: a[j] < 0}. Performing substitution alpar@9: * x[j] = 1 - x~[j] for all j in Jn, we have: alpar@9: * alpar@9: * sum a[j] x[j] >= b ==> alpar@9: * j in J alpar@9: * alpar@9: * sum a[j] x[j] + sum a[j] x[j] >= b ==> alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * sum a[j] x[j] + sum a[j] (1 - x~[j]) >= b ==> alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * sum m a[j] x[j] - sum a[j] x~[j] >= b - sum a[j]. alpar@9: * j in Jp j in Jn j in Jn alpar@9: * alpar@9: * Thus, meaning the transformation above, we can assume that in alpar@9: * inequality (1) all coefficients a[j] are positive. Moreover, we can alpar@9: * assume that b > 0, because otherwise the inequality (1) would be alpar@9: * redundant (see the routine npp_analyze_row). It is then obvious that alpar@9: * constraint (1) is equivalent to covering inequality only if: alpar@9: * alpar@9: * a[j] >= b, (2) alpar@9: * alpar@9: * for all j in J. alpar@9: * alpar@9: * Once the original inequality (1) is replaced by equivalent covering alpar@9: * inequality, we need to perform back substitution x~[j] = 1 - x[j] for alpar@9: * all j in Jn (see above). alpar@9: * alpar@9: * RECOVERING SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: static int hidden_covering(NPP *npp, struct elem *ptr, double *_b) alpar@9: { /* process inequality constraint: sum a[j] x[j] >= b; alpar@9: 0 - specified row is NOT hidden covering inequality; alpar@9: 1 - specified row is covering inequality; alpar@9: 2 - specified row is hidden covering inequality. */ alpar@9: struct elem *e; alpar@9: int neg; alpar@9: double b = *_b, eps; alpar@9: xassert(npp == npp); alpar@9: /* a[j] must be non-zero, x[j] must be binary, for all j in J */ alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: { xassert(e->aj != 0.0); alpar@9: xassert(e->xj->is_int); alpar@9: xassert(e->xj->lb == 0.0 && e->xj->ub == 1.0); alpar@9: } alpar@9: /* check if the specified inequality constraint already has the alpar@9: form of covering inequality */ alpar@9: neg = 0; /* neg is |Jn| */ alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: { if (e->aj == +1.0) alpar@9: ; alpar@9: else if (e->aj == -1.0) alpar@9: neg++; alpar@9: else alpar@9: break; alpar@9: } alpar@9: if (e == NULL) alpar@9: { /* all coefficients a[j] are +1 or -1; check rhs b */ alpar@9: if (b == (double)(1 - neg)) alpar@9: { /* it is covering inequality; no processing is needed */ alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: /* substitute x[j] = 1 - x~[j] for all j in Jn to make all a[j] alpar@9: positive; the result is a~[j] = |a[j]| and new rhs b */ alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: if (e->aj < 0) b -= e->aj; alpar@9: /* now a[j] > 0 for all j in J (actually |a[j]| are used) */ alpar@9: /* if b <= 0, skip processing--this case must not appear */ alpar@9: if (b < 1e-3) return 0; alpar@9: /* now a[j] > 0 for all j in J, and b > 0 */ alpar@9: /* the specified constraint is equivalent to covering inequality alpar@9: iff a[j] >= b for all j in J */ alpar@9: eps = 1e-9 + 1e-12 * fabs(b); alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: if (fabs(e->aj) < b - eps) return 0; alpar@9: /* perform back substitution x~[j] = 1 - x[j] and construct the alpar@9: final equivalent covering inequality in generalized format */ alpar@9: b = 1.0; alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: { if (e->aj > 0.0) alpar@9: e->aj = +1.0; alpar@9: else /* e->aj < 0.0 */ alpar@9: e->aj = -1.0, b -= 1.0; alpar@9: } alpar@9: *_b = b; alpar@9: return 2; alpar@9: } alpar@9: alpar@9: int npp_hidden_covering(NPP *npp, NPPROW *row) alpar@9: { /* identify hidden covering inequality */ alpar@9: NPPROW *copy; alpar@9: NPPAIJ *aij; alpar@9: struct elem *ptr, *e; alpar@9: int kase, ret, count = 0; alpar@9: double b; alpar@9: /* the row must be inequality constraint */ alpar@9: xassert(row->lb < row->ub); alpar@9: for (kase = 0; kase <= 1; kase++) alpar@9: { if (kase == 0) alpar@9: { /* process row lower bound */ alpar@9: if (row->lb == -DBL_MAX) continue; alpar@9: ptr = copy_form(npp, row, +1.0); alpar@9: b = + row->lb; alpar@9: } alpar@9: else alpar@9: { /* process row upper bound */ alpar@9: if (row->ub == +DBL_MAX) continue; alpar@9: ptr = copy_form(npp, row, -1.0); alpar@9: b = - row->ub; alpar@9: } alpar@9: /* now the inequality has the form "sum a[j] x[j] >= b" */ alpar@9: ret = hidden_covering(npp, ptr, &b); alpar@9: xassert(0 <= ret && ret <= 2); alpar@9: if (kase == 1 && ret == 1 || ret == 2) alpar@9: { /* the original inequality has been identified as hidden alpar@9: covering inequality */ alpar@9: count++; alpar@9: #ifdef GLP_DEBUG alpar@9: xprintf("Original constraint:\n"); alpar@9: for (aij = row->ptr; aij != NULL; aij = aij->r_next) alpar@9: xprintf(" %+g x%d", aij->val, aij->col->j); alpar@9: if (row->lb != -DBL_MAX) xprintf(", >= %g", row->lb); alpar@9: if (row->ub != +DBL_MAX) xprintf(", <= %g", row->ub); alpar@9: xprintf("\n"); alpar@9: xprintf("Equivalent covering inequality:\n"); alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: xprintf(" %sx%d", e->aj > 0.0 ? "+" : "-", e->xj->j); alpar@9: xprintf(", >= %g\n", b); alpar@9: #endif alpar@9: if (row->lb == -DBL_MAX || row->ub == +DBL_MAX) alpar@9: { /* the original row is single-sided inequality; no copy alpar@9: is needed */ alpar@9: copy = NULL; alpar@9: } alpar@9: else alpar@9: { /* the original row is double-sided inequality; we need alpar@9: to create its copy for other bound before replacing it alpar@9: with the equivalent inequality */ alpar@9: copy = npp_add_row(npp); alpar@9: if (kase == 0) alpar@9: { /* the copy is for upper bound */ alpar@9: copy->lb = -DBL_MAX, copy->ub = row->ub; alpar@9: } alpar@9: else alpar@9: { /* the copy is for lower bound */ alpar@9: copy->lb = row->lb, copy->ub = +DBL_MAX; alpar@9: } alpar@9: /* copy original row coefficients */ alpar@9: for (aij = row->ptr; aij != NULL; aij = aij->r_next) alpar@9: npp_add_aij(npp, copy, aij->col, aij->val); alpar@9: } alpar@9: /* replace the original inequality by equivalent one */ alpar@9: npp_erase_row(npp, row); alpar@9: row->lb = b, row->ub = +DBL_MAX; alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: npp_add_aij(npp, row, e->xj, e->aj); alpar@9: /* continue processing upper bound for the copy */ alpar@9: if (copy != NULL) row = copy; alpar@9: } alpar@9: drop_form(npp, ptr); alpar@9: } alpar@9: return count; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_is_partitioning - test if constraint is partitioning equality alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_is_partitioning(NPP *npp, NPPROW *row); alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * If the specified row (constraint) is partitioning equality (see alpar@9: * below), the routine npp_is_partitioning returns non-zero. Otherwise, alpar@9: * it returns zero. alpar@9: * alpar@9: * PARTITIONING EQUALITIES alpar@9: * alpar@9: * In canonical format the partitioning equality is the following: alpar@9: * alpar@9: * sum x[j] = 1, (1) alpar@9: * j in J alpar@9: * alpar@9: * where all variables x[j] are binary. This equality expresses the alpar@9: * condition that in any integer feasible solution exactly one variable alpar@9: * in set J must take non-zero (unity) value while other variables must alpar@9: * be equal to zero. W.l.o.g. it is assumed that |J| >= 2, because if alpar@9: * J is empty, the inequality (1) is infeasible, and if |J| = 1, the alpar@9: * inequality (1) is a fixing row. alpar@9: * alpar@9: * In general case the partitioning equality may include original alpar@9: * variables x[j] as well as their complements x~[j]: alpar@9: * alpar@9: * sum x[j] + sum x~[j] = 1, (2) alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * where Jp and Jn are not intersected. Therefore, using substitution alpar@9: * x~[j] = 1 - x[j] leads to the partitioning equality in generalized alpar@9: * format: alpar@9: * alpar@9: * sum x[j] - sum x[j] = 1 - |Jn|. (3) alpar@9: * j in Jp j in Jn */ alpar@9: alpar@9: int npp_is_partitioning(NPP *npp, NPPROW *row) alpar@9: { /* test if constraint is partitioning equality */ alpar@9: NPPCOL *col; alpar@9: NPPAIJ *aij; alpar@9: int b; alpar@9: xassert(npp == npp); alpar@9: if (row->lb != row->ub) return 0; alpar@9: b = 1; alpar@9: for (aij = row->ptr; aij != NULL; aij = aij->r_next) alpar@9: { col = aij->col; alpar@9: if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0)) alpar@9: return 0; alpar@9: if (aij->val == +1.0) alpar@9: ; alpar@9: else if (aij->val == -1.0) alpar@9: b--; alpar@9: else alpar@9: return 0; alpar@9: } alpar@9: if (row->lb != (double)b) return 0; alpar@9: return 1; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_reduce_ineq_coef - reduce inequality constraint coefficients alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_reduce_ineq_coef(NPP *npp, NPPROW *row); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_reduce_ineq_coef processes specified inequality alpar@9: * constraint attempting to replace it by an equivalent constraint, alpar@9: * where magnitude of coefficients at binary variables is smaller than alpar@9: * in the original constraint. If the inequality is double-sided, it is alpar@9: * replaced by a pair of single-sided inequalities, if necessary. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * The routine npp_reduce_ineq_coef returns the number of coefficients alpar@9: * reduced. alpar@9: * alpar@9: * BACKGROUND alpar@9: * alpar@9: * Consider an inequality constraint: alpar@9: * alpar@9: * sum a[j] x[j] >= b. (1) alpar@9: * j in J alpar@9: * alpar@9: * (In case of '<=' inequality it can be transformed to '>=' format by alpar@9: * multiplying both its sides by -1.) Let x[k] be a binary variable; alpar@9: * other variables can be integer as well as continuous. We can write alpar@9: * constraint (1) as follows: alpar@9: * alpar@9: * a[k] x[k] + t[k] >= b, (2) alpar@9: * alpar@9: * where: alpar@9: * alpar@9: * t[k] = sum a[j] x[j]. (3) alpar@9: * j in J\{k} alpar@9: * alpar@9: * Since x[k] is binary, constraint (2) is equivalent to disjunction of alpar@9: * the following two constraints: alpar@9: * alpar@9: * x[k] = 0, t[k] >= b (4) alpar@9: * alpar@9: * OR alpar@9: * alpar@9: * x[k] = 1, t[k] >= b - a[k]. (5) alpar@9: * alpar@9: * Let also that for the partial sum t[k] be known some its implied alpar@9: * lower bound inf t[k]. alpar@9: * alpar@9: * Case a[k] > 0. Let inf t[k] < b, since otherwise both constraints alpar@9: * (4) and (5) and therefore constraint (2) are redundant. alpar@9: * If inf t[k] > b - a[k], only constraint (5) is redundant, in which alpar@9: * case it can be replaced with the following redundant and therefore alpar@9: * equivalent constraint: alpar@9: * alpar@9: * t[k] >= b - a'[k] = inf t[k], (6) alpar@9: * alpar@9: * where: alpar@9: * alpar@9: * a'[k] = b - inf t[k]. (7) alpar@9: * alpar@9: * Thus, the original constraint (2) is equivalent to the following alpar@9: * constraint with coefficient at variable x[k] changed: alpar@9: * alpar@9: * a'[k] x[k] + t[k] >= b. (8) alpar@9: * alpar@9: * From inf t[k] < b it follows that a'[k] > 0, i.e. the coefficient alpar@9: * at x[k] keeps its sign. And from inf t[k] > b - a[k] it follows that alpar@9: * a'[k] < a[k], i.e. the coefficient reduces in magnitude. alpar@9: * alpar@9: * Case a[k] < 0. Let inf t[k] < b - a[k], since otherwise both alpar@9: * constraints (4) and (5) and therefore constraint (2) are redundant. alpar@9: * If inf t[k] > b, only constraint (4) is redundant, in which case it alpar@9: * can be replaced with the following redundant and therefore equivalent alpar@9: * constraint: alpar@9: * alpar@9: * t[k] >= b' = inf t[k]. (9) alpar@9: * alpar@9: * Rewriting constraint (5) as follows: alpar@9: * alpar@9: * t[k] >= b - a[k] = b' - a'[k], (10) alpar@9: * alpar@9: * where: alpar@9: * alpar@9: * a'[k] = a[k] + b' - b = a[k] + inf t[k] - b, (11) alpar@9: * alpar@9: * we can see that disjunction of constraint (9) and (10) is equivalent alpar@9: * to disjunction of constraint (4) and (5), from which it follows that alpar@9: * the original constraint (2) is equivalent to the following constraint alpar@9: * with both coefficient at variable x[k] and right-hand side changed: alpar@9: * alpar@9: * a'[k] x[k] + t[k] >= b'. (12) alpar@9: * alpar@9: * From inf t[k] < b - a[k] it follows that a'[k] < 0, i.e. the alpar@9: * coefficient at x[k] keeps its sign. And from inf t[k] > b it follows alpar@9: * that a'[k] > a[k], i.e. the coefficient reduces in magnitude. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * In the routine npp_reduce_ineq_coef the following implied lower alpar@9: * bound of the partial sum (3) is used: alpar@9: * alpar@9: * inf t[k] = sum a[j] l[j] + sum a[j] u[j], (13) alpar@9: * j in Jp\{k} k in Jn\{k} alpar@9: * alpar@9: * where Jp = {j : a[j] > 0}, Jn = {j : a[j] < 0}, l[j] and u[j] are alpar@9: * lower and upper bounds, resp., of variable x[j]. alpar@9: * alpar@9: * In order to compute inf t[k] more efficiently, the following formula, alpar@9: * which is equivalent to (13), is actually used: alpar@9: * alpar@9: * ( h - a[k] l[k] = h, if a[k] > 0, alpar@9: * inf t[k] = < (14) alpar@9: * ( h - a[k] u[k] = h - a[k], if a[k] < 0, alpar@9: * alpar@9: * where: alpar@9: * alpar@9: * h = sum a[j] l[j] + sum a[j] u[j] (15) alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * is the implied lower bound of row (1). alpar@9: * alpar@9: * Reduction of positive coefficient (a[k] > 0) does not change value alpar@9: * of h, since l[k] = 0. In case of reduction of negative coefficient alpar@9: * (a[k] < 0) from (11) it follows that: alpar@9: * alpar@9: * delta a[k] = a'[k] - a[k] = inf t[k] - b (> 0), (16) alpar@9: * alpar@9: * so new value of h (accounting that u[k] = 1) can be computed as alpar@9: * follows: alpar@9: * alpar@9: * h := h + delta a[k] = h + (inf t[k] - b). (17) alpar@9: * alpar@9: * RECOVERING SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: static int reduce_ineq_coef(NPP *npp, struct elem *ptr, double *_b) alpar@9: { /* process inequality constraint: sum a[j] x[j] >= b */ alpar@9: /* returns: the number of coefficients reduced */ alpar@9: struct elem *e; alpar@9: int count = 0; alpar@9: double h, inf_t, new_a, b = *_b; alpar@9: xassert(npp == npp); alpar@9: /* compute h; see (15) */ alpar@9: h = 0.0; alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: { if (e->aj > 0.0) alpar@9: { if (e->xj->lb == -DBL_MAX) goto done; alpar@9: h += e->aj * e->xj->lb; alpar@9: } alpar@9: else /* e->aj < 0.0 */ alpar@9: { if (e->xj->ub == +DBL_MAX) goto done; alpar@9: h += e->aj * e->xj->ub; alpar@9: } alpar@9: } alpar@9: /* perform reduction of coefficients at binary variables */ alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: { /* skip non-binary variable */ alpar@9: if (!(e->xj->is_int && e->xj->lb == 0.0 && e->xj->ub == 1.0)) alpar@9: continue; alpar@9: if (e->aj > 0.0) alpar@9: { /* compute inf t[k]; see (14) */ alpar@9: inf_t = h; alpar@9: if (b - e->aj < inf_t && inf_t < b) alpar@9: { /* compute reduced coefficient a'[k]; see (7) */ alpar@9: new_a = b - inf_t; alpar@9: if (new_a >= +1e-3 && alpar@9: e->aj - new_a >= 0.01 * (1.0 + e->aj)) alpar@9: { /* accept a'[k] */ alpar@9: #ifdef GLP_DEBUG alpar@9: xprintf("+"); alpar@9: #endif alpar@9: e->aj = new_a; alpar@9: count++; alpar@9: } alpar@9: } alpar@9: } alpar@9: else /* e->aj < 0.0 */ alpar@9: { /* compute inf t[k]; see (14) */ alpar@9: inf_t = h - e->aj; alpar@9: if (b < inf_t && inf_t < b - e->aj) alpar@9: { /* compute reduced coefficient a'[k]; see (11) */ alpar@9: new_a = e->aj + (inf_t - b); alpar@9: if (new_a <= -1e-3 && alpar@9: new_a - e->aj >= 0.01 * (1.0 - e->aj)) alpar@9: { /* accept a'[k] */ alpar@9: #ifdef GLP_DEBUG alpar@9: xprintf("-"); alpar@9: #endif alpar@9: e->aj = new_a; alpar@9: /* update h; see (17) */ alpar@9: h += (inf_t - b); alpar@9: /* compute b'; see (9) */ alpar@9: b = inf_t; alpar@9: count++; alpar@9: } alpar@9: } alpar@9: } alpar@9: } alpar@9: *_b = b; alpar@9: done: return count; alpar@9: } alpar@9: alpar@9: int npp_reduce_ineq_coef(NPP *npp, NPPROW *row) alpar@9: { /* reduce inequality constraint coefficients */ alpar@9: NPPROW *copy; alpar@9: NPPAIJ *aij; alpar@9: struct elem *ptr, *e; alpar@9: int kase, count[2]; alpar@9: double b; alpar@9: /* the row must be inequality constraint */ alpar@9: xassert(row->lb < row->ub); alpar@9: count[0] = count[1] = 0; alpar@9: for (kase = 0; kase <= 1; kase++) alpar@9: { if (kase == 0) alpar@9: { /* process row lower bound */ alpar@9: if (row->lb == -DBL_MAX) continue; alpar@9: #ifdef GLP_DEBUG alpar@9: xprintf("L"); alpar@9: #endif alpar@9: ptr = copy_form(npp, row, +1.0); alpar@9: b = + row->lb; alpar@9: } alpar@9: else alpar@9: { /* process row upper bound */ alpar@9: if (row->ub == +DBL_MAX) continue; alpar@9: #ifdef GLP_DEBUG alpar@9: xprintf("U"); alpar@9: #endif alpar@9: ptr = copy_form(npp, row, -1.0); alpar@9: b = - row->ub; alpar@9: } alpar@9: /* now the inequality has the form "sum a[j] x[j] >= b" */ alpar@9: count[kase] = reduce_ineq_coef(npp, ptr, &b); alpar@9: if (count[kase] > 0) alpar@9: { /* the original inequality has been replaced by equivalent alpar@9: one with coefficients reduced */ alpar@9: if (row->lb == -DBL_MAX || row->ub == +DBL_MAX) alpar@9: { /* the original row is single-sided inequality; no copy alpar@9: is needed */ alpar@9: copy = NULL; alpar@9: } alpar@9: else alpar@9: { /* the original row is double-sided inequality; we need alpar@9: to create its copy for other bound before replacing it alpar@9: with the equivalent inequality */ alpar@9: #ifdef GLP_DEBUG alpar@9: xprintf("*"); alpar@9: #endif alpar@9: copy = npp_add_row(npp); alpar@9: if (kase == 0) alpar@9: { /* the copy is for upper bound */ alpar@9: copy->lb = -DBL_MAX, copy->ub = row->ub; alpar@9: } alpar@9: else alpar@9: { /* the copy is for lower bound */ alpar@9: copy->lb = row->lb, copy->ub = +DBL_MAX; alpar@9: } alpar@9: /* copy original row coefficients */ alpar@9: for (aij = row->ptr; aij != NULL; aij = aij->r_next) alpar@9: npp_add_aij(npp, copy, aij->col, aij->val); alpar@9: } alpar@9: /* replace the original inequality by equivalent one */ alpar@9: npp_erase_row(npp, row); alpar@9: row->lb = b, row->ub = +DBL_MAX; alpar@9: for (e = ptr; e != NULL; e = e->next) alpar@9: npp_add_aij(npp, row, e->xj, e->aj); alpar@9: /* continue processing upper bound for the copy */ alpar@9: if (copy != NULL) row = copy; alpar@9: } alpar@9: drop_form(npp, ptr); alpar@9: } alpar@9: return count[0] + count[1]; alpar@9: } alpar@9: alpar@9: /* eof */