lemon-project-template-glpk
diff deps/glpk/src/glpnpp06.c @ 9:33de93886c88
Import GLPK 4.47
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 20:59:10 +0100 |
parents | |
children |
line diff
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/deps/glpk/src/glpnpp06.c Sun Nov 06 20:59:10 2011 +0100 1.3 @@ -0,0 +1,1500 @@ 1.4 +/* glpnpp06.c (translate feasibility problem to CNF-SAT) */ 1.5 + 1.6 +/*********************************************************************** 1.7 +* This code is part of GLPK (GNU Linear Programming Kit). 1.8 +* 1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 1.10 +* 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, 1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved. 1.12 +* E-mail: <mao@gnu.org>. 1.13 +* 1.14 +* GLPK is free software: you can redistribute it and/or modify it 1.15 +* under the terms of the GNU General Public License as published by 1.16 +* the Free Software Foundation, either version 3 of the License, or 1.17 +* (at your option) any later version. 1.18 +* 1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT 1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public 1.22 +* License for more details. 1.23 +* 1.24 +* You should have received a copy of the GNU General Public License 1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>. 1.26 +***********************************************************************/ 1.27 + 1.28 +#include "glpnpp.h" 1.29 + 1.30 +/*********************************************************************** 1.31 +* npp_sat_free_row - process free (unbounded) row 1.32 +* 1.33 +* This routine processes row p, which is free (i.e. has no finite 1.34 +* bounds): 1.35 +* 1.36 +* -inf < sum a[p,j] x[j] < +inf. (1) 1.37 +* 1.38 +* The constraint (1) cannot be active and therefore it is redundant, 1.39 +* so the routine simply removes it from the original problem. */ 1.40 + 1.41 +void npp_sat_free_row(NPP *npp, NPPROW *p) 1.42 +{ /* the row should be free */ 1.43 + xassert(p->lb == -DBL_MAX && p->ub == +DBL_MAX); 1.44 + /* remove the row from the problem */ 1.45 + npp_del_row(npp, p); 1.46 + return; 1.47 +} 1.48 + 1.49 +/*********************************************************************** 1.50 +* npp_sat_fixed_col - process fixed column 1.51 +* 1.52 +* This routine processes column q, which is fixed: 1.53 +* 1.54 +* x[q] = s[q], (1) 1.55 +* 1.56 +* where s[q] is a fixed column value. 1.57 +* 1.58 +* The routine substitutes fixed value s[q] into constraint rows and 1.59 +* then removes column x[q] from the original problem. 1.60 +* 1.61 +* Substitution of x[q] = s[q] into row i gives: 1.62 +* 1.63 +* L[i] <= sum a[i,j] x[j] <= U[i] ==> 1.64 +* j 1.65 +* 1.66 +* L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==> 1.67 +* j!=q 1.68 +* 1.69 +* L[i] <= sum a[i,j] x[j] + a[i,q] s[q] <= U[i] ==> 1.70 +* j!=q 1.71 +* 1.72 +* L~[i] <= sum a[i,j] x[j] <= U~[i], 1.73 +* j!=q 1.74 +* 1.75 +* where 1.76 +* 1.77 +* L~[i] = L[i] - a[i,q] s[q], (2) 1.78 +* 1.79 +* U~[i] = U[i] - a[i,q] s[q] (3) 1.80 +* 1.81 +* are, respectively, lower and upper bound of row i in the transformed 1.82 +* problem. 1.83 +* 1.84 +* On recovering solution x[q] is assigned the value of s[q]. */ 1.85 + 1.86 +struct sat_fixed_col 1.87 +{ /* fixed column */ 1.88 + int q; 1.89 + /* column reference number for variable x[q] */ 1.90 + int s; 1.91 + /* value, at which x[q] is fixed */ 1.92 +}; 1.93 + 1.94 +static int rcv_sat_fixed_col(NPP *, void *); 1.95 + 1.96 +int npp_sat_fixed_col(NPP *npp, NPPCOL *q) 1.97 +{ struct sat_fixed_col *info; 1.98 + NPPROW *i; 1.99 + NPPAIJ *aij; 1.100 + int temp; 1.101 + /* the column should be fixed */ 1.102 + xassert(q->lb == q->ub); 1.103 + /* create transformation stack entry */ 1.104 + info = npp_push_tse(npp, 1.105 + rcv_sat_fixed_col, sizeof(struct sat_fixed_col)); 1.106 + info->q = q->j; 1.107 + info->s = (int)q->lb; 1.108 + xassert((double)info->s == q->lb); 1.109 + /* substitute x[q] = s[q] into constraint rows */ 1.110 + if (info->s == 0) 1.111 + goto skip; 1.112 + for (aij = q->ptr; aij != NULL; aij = aij->c_next) 1.113 + { i = aij->row; 1.114 + if (i->lb != -DBL_MAX) 1.115 + { i->lb -= aij->val * (double)info->s; 1.116 + temp = (int)i->lb; 1.117 + if ((double)temp != i->lb) 1.118 + return 1; /* integer arithmetic error */ 1.119 + } 1.120 + if (i->ub != +DBL_MAX) 1.121 + { i->ub -= aij->val * (double)info->s; 1.122 + temp = (int)i->ub; 1.123 + if ((double)temp != i->ub) 1.124 + return 2; /* integer arithmetic error */ 1.125 + } 1.126 + } 1.127 +skip: /* remove the column from the problem */ 1.128 + npp_del_col(npp, q); 1.129 + return 0; 1.130 +} 1.131 + 1.132 +static int rcv_sat_fixed_col(NPP *npp, void *info_) 1.133 +{ struct sat_fixed_col *info = info_; 1.134 + npp->c_value[info->q] = (double)info->s; 1.135 + return 0; 1.136 +} 1.137 + 1.138 +/*********************************************************************** 1.139 +* npp_sat_is_bin_comb - test if row is binary combination 1.140 +* 1.141 +* This routine tests if the specified row is a binary combination, 1.142 +* i.e. all its constraint coefficients are +1 and -1 and all variables 1.143 +* are binary. If the test was passed, the routine returns non-zero, 1.144 +* otherwise zero. */ 1.145 + 1.146 +int npp_sat_is_bin_comb(NPP *npp, NPPROW *row) 1.147 +{ NPPCOL *col; 1.148 + NPPAIJ *aij; 1.149 + xassert(npp == npp); 1.150 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.151 + { if (!(aij->val == +1.0 || aij->val == -1.0)) 1.152 + return 0; /* non-unity coefficient */ 1.153 + col = aij->col; 1.154 + if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0)) 1.155 + return 0; /* non-binary column */ 1.156 + } 1.157 + return 1; /* test was passed */ 1.158 +} 1.159 + 1.160 +/*********************************************************************** 1.161 +* npp_sat_num_pos_coef - determine number of positive coefficients 1.162 +* 1.163 +* This routine returns the number of positive coefficients in the 1.164 +* specified row. */ 1.165 + 1.166 +int npp_sat_num_pos_coef(NPP *npp, NPPROW *row) 1.167 +{ NPPAIJ *aij; 1.168 + int num = 0; 1.169 + xassert(npp == npp); 1.170 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.171 + { if (aij->val > 0.0) 1.172 + num++; 1.173 + } 1.174 + return num; 1.175 +} 1.176 + 1.177 +/*********************************************************************** 1.178 +* npp_sat_num_neg_coef - determine number of negative coefficients 1.179 +* 1.180 +* This routine returns the number of negative coefficients in the 1.181 +* specified row. */ 1.182 + 1.183 +int npp_sat_num_neg_coef(NPP *npp, NPPROW *row) 1.184 +{ NPPAIJ *aij; 1.185 + int num = 0; 1.186 + xassert(npp == npp); 1.187 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.188 + { if (aij->val < 0.0) 1.189 + num++; 1.190 + } 1.191 + return num; 1.192 +} 1.193 + 1.194 +/*********************************************************************** 1.195 +* npp_sat_is_cover_ineq - test if row is covering inequality 1.196 +* 1.197 +* The canonical form of a covering inequality is the following: 1.198 +* 1.199 +* sum x[j] >= 1, (1) 1.200 +* j in J 1.201 +* 1.202 +* where all x[j] are binary variables. 1.203 +* 1.204 +* In general case a covering inequality may have one of the following 1.205 +* two forms: 1.206 +* 1.207 +* sum x[j] - sum x[j] >= 1 - |J-|, (2) 1.208 +* j in J+ j in J- 1.209 +* 1.210 +* 1.211 +* sum x[j] - sum x[j] <= |J+| - 1. (3) 1.212 +* j in J+ j in J- 1.213 +* 1.214 +* Obviously, the inequality (2) can be transformed to the form (1) by 1.215 +* substitution x[j] = 1 - x'[j] for all j in J-, where x'[j] is the 1.216 +* negation of variable x[j]. And the inequality (3) can be transformed 1.217 +* to (2) by multiplying both left- and right-hand sides by -1. 1.218 +* 1.219 +* This routine returns one of the following codes: 1.220 +* 1.221 +* 0, if the specified row is not a covering inequality; 1.222 +* 1.223 +* 1, if the specified row has the form (2); 1.224 +* 1.225 +* 2, if the specified row has the form (3). */ 1.226 + 1.227 +int npp_sat_is_cover_ineq(NPP *npp, NPPROW *row) 1.228 +{ xassert(npp == npp); 1.229 + if (row->lb != -DBL_MAX && row->ub == +DBL_MAX) 1.230 + { /* row is inequality of '>=' type */ 1.231 + if (npp_sat_is_bin_comb(npp, row)) 1.232 + { /* row is a binary combination */ 1.233 + if (row->lb == 1.0 - npp_sat_num_neg_coef(npp, row)) 1.234 + { /* row has the form (2) */ 1.235 + return 1; 1.236 + } 1.237 + } 1.238 + } 1.239 + else if (row->lb == -DBL_MAX && row->ub != +DBL_MAX) 1.240 + { /* row is inequality of '<=' type */ 1.241 + if (npp_sat_is_bin_comb(npp, row)) 1.242 + { /* row is a binary combination */ 1.243 + if (row->ub == npp_sat_num_pos_coef(npp, row) - 1.0) 1.244 + { /* row has the form (3) */ 1.245 + return 2; 1.246 + } 1.247 + } 1.248 + } 1.249 + /* row is not a covering inequality */ 1.250 + return 0; 1.251 +} 1.252 + 1.253 +/*********************************************************************** 1.254 +* npp_sat_is_pack_ineq - test if row is packing inequality 1.255 +* 1.256 +* The canonical form of a packing inequality is the following: 1.257 +* 1.258 +* sum x[j] <= 1, (1) 1.259 +* j in J 1.260 +* 1.261 +* where all x[j] are binary variables. 1.262 +* 1.263 +* In general case a packing inequality may have one of the following 1.264 +* two forms: 1.265 +* 1.266 +* sum x[j] - sum x[j] <= 1 - |J-|, (2) 1.267 +* j in J+ j in J- 1.268 +* 1.269 +* 1.270 +* sum x[j] - sum x[j] >= |J+| - 1. (3) 1.271 +* j in J+ j in J- 1.272 +* 1.273 +* Obviously, the inequality (2) can be transformed to the form (1) by 1.274 +* substitution x[j] = 1 - x'[j] for all j in J-, where x'[j] is the 1.275 +* negation of variable x[j]. And the inequality (3) can be transformed 1.276 +* to (2) by multiplying both left- and right-hand sides by -1. 1.277 +* 1.278 +* This routine returns one of the following codes: 1.279 +* 1.280 +* 0, if the specified row is not a packing inequality; 1.281 +* 1.282 +* 1, if the specified row has the form (2); 1.283 +* 1.284 +* 2, if the specified row has the form (3). */ 1.285 + 1.286 +int npp_sat_is_pack_ineq(NPP *npp, NPPROW *row) 1.287 +{ xassert(npp == npp); 1.288 + if (row->lb == -DBL_MAX && row->ub != +DBL_MAX) 1.289 + { /* row is inequality of '<=' type */ 1.290 + if (npp_sat_is_bin_comb(npp, row)) 1.291 + { /* row is a binary combination */ 1.292 + if (row->ub == 1.0 - npp_sat_num_neg_coef(npp, row)) 1.293 + { /* row has the form (2) */ 1.294 + return 1; 1.295 + } 1.296 + } 1.297 + } 1.298 + else if (row->lb != -DBL_MAX && row->ub == +DBL_MAX) 1.299 + { /* row is inequality of '>=' type */ 1.300 + if (npp_sat_is_bin_comb(npp, row)) 1.301 + { /* row is a binary combination */ 1.302 + if (row->lb == npp_sat_num_pos_coef(npp, row) - 1.0) 1.303 + { /* row has the form (3) */ 1.304 + return 2; 1.305 + } 1.306 + } 1.307 + } 1.308 + /* row is not a packing inequality */ 1.309 + return 0; 1.310 +} 1.311 + 1.312 +/*********************************************************************** 1.313 +* npp_sat_is_partn_eq - test if row is partitioning equality 1.314 +* 1.315 +* The canonical form of a partitioning equality is the following: 1.316 +* 1.317 +* sum x[j] = 1, (1) 1.318 +* j in J 1.319 +* 1.320 +* where all x[j] are binary variables. 1.321 +* 1.322 +* In general case a partitioning equality may have one of the following 1.323 +* two forms: 1.324 +* 1.325 +* sum x[j] - sum x[j] = 1 - |J-|, (2) 1.326 +* j in J+ j in J- 1.327 +* 1.328 +* 1.329 +* sum x[j] - sum x[j] = |J+| - 1. (3) 1.330 +* j in J+ j in J- 1.331 +* 1.332 +* Obviously, the equality (2) can be transformed to the form (1) by 1.333 +* substitution x[j] = 1 - x'[j] for all j in J-, where x'[j] is the 1.334 +* negation of variable x[j]. And the equality (3) can be transformed 1.335 +* to (2) by multiplying both left- and right-hand sides by -1. 1.336 +* 1.337 +* This routine returns one of the following codes: 1.338 +* 1.339 +* 0, if the specified row is not a partitioning equality; 1.340 +* 1.341 +* 1, if the specified row has the form (2); 1.342 +* 1.343 +* 2, if the specified row has the form (3). */ 1.344 + 1.345 +int npp_sat_is_partn_eq(NPP *npp, NPPROW *row) 1.346 +{ xassert(npp == npp); 1.347 + if (row->lb == row->ub) 1.348 + { /* row is equality constraint */ 1.349 + if (npp_sat_is_bin_comb(npp, row)) 1.350 + { /* row is a binary combination */ 1.351 + if (row->lb == 1.0 - npp_sat_num_neg_coef(npp, row)) 1.352 + { /* row has the form (2) */ 1.353 + return 1; 1.354 + } 1.355 + if (row->ub == npp_sat_num_pos_coef(npp, row) - 1.0) 1.356 + { /* row has the form (3) */ 1.357 + return 2; 1.358 + } 1.359 + } 1.360 + } 1.361 + /* row is not a partitioning equality */ 1.362 + return 0; 1.363 +} 1.364 + 1.365 +/*********************************************************************** 1.366 +* npp_sat_reverse_row - multiply both sides of row by -1 1.367 +* 1.368 +* This routines multiplies by -1 both left- and right-hand sides of 1.369 +* the specified row: 1.370 +* 1.371 +* L <= sum x[j] <= U, 1.372 +* 1.373 +* that results in the following row: 1.374 +* 1.375 +* -U <= sum (-x[j]) <= -L. 1.376 +* 1.377 +* If no integer overflow occured, the routine returns zero, otherwise 1.378 +* non-zero. */ 1.379 + 1.380 +int npp_sat_reverse_row(NPP *npp, NPPROW *row) 1.381 +{ NPPAIJ *aij; 1.382 + int temp, ret = 0; 1.383 + double old_lb, old_ub; 1.384 + xassert(npp == npp); 1.385 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.386 + { aij->val = -aij->val; 1.387 + temp = (int)aij->val; 1.388 + if ((double)temp != aij->val) 1.389 + ret = 1; 1.390 + } 1.391 + old_lb = row->lb, old_ub = row->ub; 1.392 + if (old_ub == +DBL_MAX) 1.393 + row->lb = -DBL_MAX; 1.394 + else 1.395 + { row->lb = -old_ub; 1.396 + temp = (int)row->lb; 1.397 + if ((double)temp != row->lb) 1.398 + ret = 2; 1.399 + } 1.400 + if (old_lb == -DBL_MAX) 1.401 + row->ub = +DBL_MAX; 1.402 + else 1.403 + { row->ub = -old_lb; 1.404 + temp = (int)row->ub; 1.405 + if ((double)temp != row->ub) 1.406 + ret = 3; 1.407 + } 1.408 + return ret; 1.409 +} 1.410 + 1.411 +/*********************************************************************** 1.412 +* npp_sat_split_pack - split packing inequality 1.413 +* 1.414 +* Let there be given a packing inequality in canonical form: 1.415 +* 1.416 +* sum t[j] <= 1, (1) 1.417 +* j in J 1.418 +* 1.419 +* where t[j] = x[j] or t[j] = 1 - x[j], x[j] is a binary variable. 1.420 +* And let J = J1 U J2 is a partition of the set of literals. Then the 1.421 +* inequality (1) is obviously equivalent to the following two packing 1.422 +* inequalities: 1.423 +* 1.424 +* sum t[j] <= y <--> sum t[j] + (1 - y) <= 1, (2) 1.425 +* j in J1 j in J1 1.426 +* 1.427 +* sum t[j] <= 1 - y <--> sum t[j] + y <= 1, (3) 1.428 +* j in J2 j in J2 1.429 +* 1.430 +* where y is a new binary variable added to the transformed problem. 1.431 +* 1.432 +* Assuming that the specified row is a packing inequality (1), this 1.433 +* routine constructs the set J1 by including there first nlit literals 1.434 +* (terms) from the specified row, and the set J2 = J \ J1. Then the 1.435 +* routine creates a new row, which corresponds to inequality (2), and 1.436 +* replaces the specified row with inequality (3). */ 1.437 + 1.438 +NPPROW *npp_sat_split_pack(NPP *npp, NPPROW *row, int nlit) 1.439 +{ NPPROW *rrr; 1.440 + NPPCOL *col; 1.441 + NPPAIJ *aij; 1.442 + int k; 1.443 + /* original row should be packing inequality (1) */ 1.444 + xassert(npp_sat_is_pack_ineq(npp, row) == 1); 1.445 + /* and nlit should be less than the number of literals (terms) 1.446 + in the original row */ 1.447 + xassert(0 < nlit && nlit < npp_row_nnz(npp, row)); 1.448 + /* create new row corresponding to inequality (2) */ 1.449 + rrr = npp_add_row(npp); 1.450 + rrr->lb = -DBL_MAX, rrr->ub = 1.0; 1.451 + /* move first nlit literals (terms) from the original row to the 1.452 + new row; the original row becomes inequality (3) */ 1.453 + for (k = 1; k <= nlit; k++) 1.454 + { aij = row->ptr; 1.455 + xassert(aij != NULL); 1.456 + /* add literal to the new row */ 1.457 + npp_add_aij(npp, rrr, aij->col, aij->val); 1.458 + /* correct rhs */ 1.459 + if (aij->val < 0.0) 1.460 + rrr->ub -= 1.0, row->ub += 1.0; 1.461 + /* remove literal from the original row */ 1.462 + npp_del_aij(npp, aij); 1.463 + } 1.464 + /* create new binary variable y */ 1.465 + col = npp_add_col(npp); 1.466 + col->is_int = 1, col->lb = 0.0, col->ub = 1.0; 1.467 + /* include literal (1 - y) in the new row */ 1.468 + npp_add_aij(npp, rrr, col, -1.0); 1.469 + rrr->ub -= 1.0; 1.470 + /* include literal y in the original row */ 1.471 + npp_add_aij(npp, row, col, +1.0); 1.472 + return rrr; 1.473 +} 1.474 + 1.475 +/*********************************************************************** 1.476 +* npp_sat_encode_pack - encode packing inequality 1.477 +* 1.478 +* Given a packing inequality in canonical form: 1.479 +* 1.480 +* sum t[j] <= 1, (1) 1.481 +* j in J 1.482 +* 1.483 +* where t[j] = x[j] or t[j] = 1 - x[j], x[j] is a binary variable, 1.484 +* this routine translates it to CNF by replacing it with the following 1.485 +* equivalent set of edge packing inequalities: 1.486 +* 1.487 +* t[j] + t[k] <= 1 for all j, k in J, j != k. (2) 1.488 +* 1.489 +* Then the routine transforms each edge packing inequality (2) to 1.490 +* corresponding covering inequality (that encodes two-literal clause) 1.491 +* by multiplying both its part by -1: 1.492 +* 1.493 +* - t[j] - t[k] >= -1 <--> (1 - t[j]) + (1 - t[k]) >= 1. (3) 1.494 +* 1.495 +* On exit the routine removes the original row from the problem. */ 1.496 + 1.497 +void npp_sat_encode_pack(NPP *npp, NPPROW *row) 1.498 +{ NPPROW *rrr; 1.499 + NPPAIJ *aij, *aik; 1.500 + /* original row should be packing inequality (1) */ 1.501 + xassert(npp_sat_is_pack_ineq(npp, row) == 1); 1.502 + /* create equivalent system of covering inequalities (3) */ 1.503 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.504 + { /* due to symmetry only one of inequalities t[j] + t[k] <= 1 1.505 + and t[k] <= t[j] <= 1 can be considered */ 1.506 + for (aik = aij->r_next; aik != NULL; aik = aik->r_next) 1.507 + { /* create edge packing inequality (2) */ 1.508 + rrr = npp_add_row(npp); 1.509 + rrr->lb = -DBL_MAX, rrr->ub = 1.0; 1.510 + npp_add_aij(npp, rrr, aij->col, aij->val); 1.511 + if (aij->val < 0.0) 1.512 + rrr->ub -= 1.0; 1.513 + npp_add_aij(npp, rrr, aik->col, aik->val); 1.514 + if (aik->val < 0.0) 1.515 + rrr->ub -= 1.0; 1.516 + /* and transform it to covering inequality (3) */ 1.517 + npp_sat_reverse_row(npp, rrr); 1.518 + xassert(npp_sat_is_cover_ineq(npp, rrr) == 1); 1.519 + } 1.520 + } 1.521 + /* remove the original row from the problem */ 1.522 + npp_del_row(npp, row); 1.523 + return; 1.524 +} 1.525 + 1.526 +/*********************************************************************** 1.527 +* npp_sat_encode_sum2 - encode 2-bit summation 1.528 +* 1.529 +* Given a set containing two literals x and y this routine encodes 1.530 +* the equality 1.531 +* 1.532 +* x + y = s + 2 * c, (1) 1.533 +* 1.534 +* where 1.535 +* 1.536 +* s = (x + y) % 2 (2) 1.537 +* 1.538 +* is a binary variable modeling the low sum bit, and 1.539 +* 1.540 +* c = (x + y) / 2 (3) 1.541 +* 1.542 +* is a binary variable modeling the high (carry) sum bit. */ 1.543 + 1.544 +void npp_sat_encode_sum2(NPP *npp, NPPLSE *set, NPPSED *sed) 1.545 +{ NPPROW *row; 1.546 + int x, y, s, c; 1.547 + /* the set should contain exactly two literals */ 1.548 + xassert(set != NULL); 1.549 + xassert(set->next != NULL); 1.550 + xassert(set->next->next == NULL); 1.551 + sed->x = set->lit; 1.552 + xassert(sed->x.neg == 0 || sed->x.neg == 1); 1.553 + sed->y = set->next->lit; 1.554 + xassert(sed->y.neg == 0 || sed->y.neg == 1); 1.555 + sed->z.col = NULL, sed->z.neg = 0; 1.556 + /* perform encoding s = (x + y) % 2 */ 1.557 + sed->s = npp_add_col(npp); 1.558 + sed->s->is_int = 1, sed->s->lb = 0.0, sed->s->ub = 1.0; 1.559 + for (x = 0; x <= 1; x++) 1.560 + { for (y = 0; y <= 1; y++) 1.561 + { for (s = 0; s <= 1; s++) 1.562 + { if ((x + y) % 2 != s) 1.563 + { /* generate CNF clause to disable infeasible 1.564 + combination */ 1.565 + row = npp_add_row(npp); 1.566 + row->lb = 1.0, row->ub = +DBL_MAX; 1.567 + if (x == sed->x.neg) 1.568 + npp_add_aij(npp, row, sed->x.col, +1.0); 1.569 + else 1.570 + { npp_add_aij(npp, row, sed->x.col, -1.0); 1.571 + row->lb -= 1.0; 1.572 + } 1.573 + if (y == sed->y.neg) 1.574 + npp_add_aij(npp, row, sed->y.col, +1.0); 1.575 + else 1.576 + { npp_add_aij(npp, row, sed->y.col, -1.0); 1.577 + row->lb -= 1.0; 1.578 + } 1.579 + if (s == 0) 1.580 + npp_add_aij(npp, row, sed->s, +1.0); 1.581 + else 1.582 + { npp_add_aij(npp, row, sed->s, -1.0); 1.583 + row->lb -= 1.0; 1.584 + } 1.585 + } 1.586 + } 1.587 + } 1.588 + } 1.589 + /* perform encoding c = (x + y) / 2 */ 1.590 + sed->c = npp_add_col(npp); 1.591 + sed->c->is_int = 1, sed->c->lb = 0.0, sed->c->ub = 1.0; 1.592 + for (x = 0; x <= 1; x++) 1.593 + { for (y = 0; y <= 1; y++) 1.594 + { for (c = 0; c <= 1; c++) 1.595 + { if ((x + y) / 2 != c) 1.596 + { /* generate CNF clause to disable infeasible 1.597 + combination */ 1.598 + row = npp_add_row(npp); 1.599 + row->lb = 1.0, row->ub = +DBL_MAX; 1.600 + if (x == sed->x.neg) 1.601 + npp_add_aij(npp, row, sed->x.col, +1.0); 1.602 + else 1.603 + { npp_add_aij(npp, row, sed->x.col, -1.0); 1.604 + row->lb -= 1.0; 1.605 + } 1.606 + if (y == sed->y.neg) 1.607 + npp_add_aij(npp, row, sed->y.col, +1.0); 1.608 + else 1.609 + { npp_add_aij(npp, row, sed->y.col, -1.0); 1.610 + row->lb -= 1.0; 1.611 + } 1.612 + if (c == 0) 1.613 + npp_add_aij(npp, row, sed->c, +1.0); 1.614 + else 1.615 + { npp_add_aij(npp, row, sed->c, -1.0); 1.616 + row->lb -= 1.0; 1.617 + } 1.618 + } 1.619 + } 1.620 + } 1.621 + } 1.622 + return; 1.623 +} 1.624 + 1.625 +/*********************************************************************** 1.626 +* npp_sat_encode_sum3 - encode 3-bit summation 1.627 +* 1.628 +* Given a set containing at least three literals this routine chooses 1.629 +* some literals x, y, z from that set and encodes the equality 1.630 +* 1.631 +* x + y + z = s + 2 * c, (1) 1.632 +* 1.633 +* where 1.634 +* 1.635 +* s = (x + y + z) % 2 (2) 1.636 +* 1.637 +* is a binary variable modeling the low sum bit, and 1.638 +* 1.639 +* c = (x + y + z) / 2 (3) 1.640 +* 1.641 +* is a binary variable modeling the high (carry) sum bit. */ 1.642 + 1.643 +void npp_sat_encode_sum3(NPP *npp, NPPLSE *set, NPPSED *sed) 1.644 +{ NPPROW *row; 1.645 + int x, y, z, s, c; 1.646 + /* the set should contain at least three literals */ 1.647 + xassert(set != NULL); 1.648 + xassert(set->next != NULL); 1.649 + xassert(set->next->next != NULL); 1.650 + sed->x = set->lit; 1.651 + xassert(sed->x.neg == 0 || sed->x.neg == 1); 1.652 + sed->y = set->next->lit; 1.653 + xassert(sed->y.neg == 0 || sed->y.neg == 1); 1.654 + sed->z = set->next->next->lit; 1.655 + xassert(sed->z.neg == 0 || sed->z.neg == 1); 1.656 + /* perform encoding s = (x + y + z) % 2 */ 1.657 + sed->s = npp_add_col(npp); 1.658 + sed->s->is_int = 1, sed->s->lb = 0.0, sed->s->ub = 1.0; 1.659 + for (x = 0; x <= 1; x++) 1.660 + { for (y = 0; y <= 1; y++) 1.661 + { for (z = 0; z <= 1; z++) 1.662 + { for (s = 0; s <= 1; s++) 1.663 + { if ((x + y + z) % 2 != s) 1.664 + { /* generate CNF clause to disable infeasible 1.665 + combination */ 1.666 + row = npp_add_row(npp); 1.667 + row->lb = 1.0, row->ub = +DBL_MAX; 1.668 + if (x == sed->x.neg) 1.669 + npp_add_aij(npp, row, sed->x.col, +1.0); 1.670 + else 1.671 + { npp_add_aij(npp, row, sed->x.col, -1.0); 1.672 + row->lb -= 1.0; 1.673 + } 1.674 + if (y == sed->y.neg) 1.675 + npp_add_aij(npp, row, sed->y.col, +1.0); 1.676 + else 1.677 + { npp_add_aij(npp, row, sed->y.col, -1.0); 1.678 + row->lb -= 1.0; 1.679 + } 1.680 + if (z == sed->z.neg) 1.681 + npp_add_aij(npp, row, sed->z.col, +1.0); 1.682 + else 1.683 + { npp_add_aij(npp, row, sed->z.col, -1.0); 1.684 + row->lb -= 1.0; 1.685 + } 1.686 + if (s == 0) 1.687 + npp_add_aij(npp, row, sed->s, +1.0); 1.688 + else 1.689 + { npp_add_aij(npp, row, sed->s, -1.0); 1.690 + row->lb -= 1.0; 1.691 + } 1.692 + } 1.693 + } 1.694 + } 1.695 + } 1.696 + } 1.697 + /* perform encoding c = (x + y + z) / 2 */ 1.698 + sed->c = npp_add_col(npp); 1.699 + sed->c->is_int = 1, sed->c->lb = 0.0, sed->c->ub = 1.0; 1.700 + for (x = 0; x <= 1; x++) 1.701 + { for (y = 0; y <= 1; y++) 1.702 + { for (z = 0; z <= 1; z++) 1.703 + { for (c = 0; c <= 1; c++) 1.704 + { if ((x + y + z) / 2 != c) 1.705 + { /* generate CNF clause to disable infeasible 1.706 + combination */ 1.707 + row = npp_add_row(npp); 1.708 + row->lb = 1.0, row->ub = +DBL_MAX; 1.709 + if (x == sed->x.neg) 1.710 + npp_add_aij(npp, row, sed->x.col, +1.0); 1.711 + else 1.712 + { npp_add_aij(npp, row, sed->x.col, -1.0); 1.713 + row->lb -= 1.0; 1.714 + } 1.715 + if (y == sed->y.neg) 1.716 + npp_add_aij(npp, row, sed->y.col, +1.0); 1.717 + else 1.718 + { npp_add_aij(npp, row, sed->y.col, -1.0); 1.719 + row->lb -= 1.0; 1.720 + } 1.721 + if (z == sed->z.neg) 1.722 + npp_add_aij(npp, row, sed->z.col, +1.0); 1.723 + else 1.724 + { npp_add_aij(npp, row, sed->z.col, -1.0); 1.725 + row->lb -= 1.0; 1.726 + } 1.727 + if (c == 0) 1.728 + npp_add_aij(npp, row, sed->c, +1.0); 1.729 + else 1.730 + { npp_add_aij(npp, row, sed->c, -1.0); 1.731 + row->lb -= 1.0; 1.732 + } 1.733 + } 1.734 + } 1.735 + } 1.736 + } 1.737 + } 1.738 + return; 1.739 +} 1.740 + 1.741 +/*********************************************************************** 1.742 +* npp_sat_encode_sum_ax - encode linear combination of 0-1 variables 1.743 +* 1.744 +* PURPOSE 1.745 +* 1.746 +* Given a linear combination of binary variables: 1.747 +* 1.748 +* sum a[j] x[j], (1) 1.749 +* j 1.750 +* 1.751 +* which is the linear form of the specified row, this routine encodes 1.752 +* (i.e. translates to CNF) the following equality: 1.753 +* 1.754 +* n 1.755 +* sum |a[j]| t[j] = sum 2**(k-1) * y[k], (2) 1.756 +* j k=1 1.757 +* 1.758 +* where t[j] = x[j] (if a[j] > 0) or t[j] = 1 - x[j] (if a[j] < 0), 1.759 +* and y[k] is either t[j] or a new literal created by the routine or 1.760 +* a constant zero. Note that the sum in the right-hand side of (2) can 1.761 +* be thought as a n-bit representation of the sum in the left-hand 1.762 +* side, which is a non-negative integer number. 1.763 +* 1.764 +* ALGORITHM 1.765 +* 1.766 +* First, the number of bits, n, sufficient to represent any value in 1.767 +* the left-hand side of (2) is determined. Obviously, n is the number 1.768 +* of bits sufficient to represent the sum (sum |a[j]|). 1.769 +* 1.770 +* Let 1.771 +* 1.772 +* n 1.773 +* |a[j]| = sum 2**(k-1) b[j,k], (3) 1.774 +* k=1 1.775 +* 1.776 +* where b[j,k] is k-th bit in a n-bit representation of |a[j]|. Then 1.777 +* 1.778 +* m n 1.779 +* sum |a[j]| * t[j] = sum 2**(k-1) sum b[j,k] * t[j]. (4) 1.780 +* j k=1 j=1 1.781 +* 1.782 +* Introducing the set 1.783 +* 1.784 +* J[k] = { j : b[j,k] = 1 } (5) 1.785 +* 1.786 +* allows rewriting (4) as follows: 1.787 +* 1.788 +* n 1.789 +* sum |a[j]| * t[j] = sum 2**(k-1) sum t[j]. (6) 1.790 +* j k=1 j in J[k] 1.791 +* 1.792 +* Thus, our goal is to provide |J[k]| <= 1 for all k, in which case 1.793 +* we will have the representation (1). 1.794 +* 1.795 +* Let |J[k]| = 2, i.e. J[k] has exactly two literals u and v. In this 1.796 +* case we can apply the following transformation: 1.797 +* 1.798 +* u + v = s + 2 * c, (7) 1.799 +* 1.800 +* where s and c are, respectively, low (sum) and high (carry) bits of 1.801 +* the sum of two bits. This allows to replace two literals u and v in 1.802 +* J[k] by one literal s, and carry out literal c to J[k+1]. 1.803 +* 1.804 +* If |J[k]| >= 3, i.e. J[k] has at least three literals u, v, and w, 1.805 +* we can apply the following transformation: 1.806 +* 1.807 +* u + v + w = s + 2 * c. (8) 1.808 +* 1.809 +* Again, literal s replaces literals u, v, and w in J[k], and literal 1.810 +* c goes into J[k+1]. 1.811 +* 1.812 +* On exit the routine stores each literal from J[k] in element y[k], 1.813 +* 1 <= k <= n. If J[k] is empty, y[k] is set to constant false. 1.814 +* 1.815 +* RETURNS 1.816 +* 1.817 +* The routine returns n, the number of literals in the right-hand side 1.818 +* of (2), 0 <= n <= NBIT_MAX. If the sum (sum |a[j]|) is too large, so 1.819 +* more than NBIT_MAX (= 31) literals are needed to encode the original 1.820 +* linear combination, the routine returns a negative value. */ 1.821 + 1.822 +#define NBIT_MAX 31 1.823 +/* maximal number of literals in the right hand-side of (2) */ 1.824 + 1.825 +static NPPLSE *remove_lse(NPP *npp, NPPLSE *set, NPPCOL *col) 1.826 +{ /* remove specified literal from specified literal set */ 1.827 + NPPLSE *lse, *prev = NULL; 1.828 + for (lse = set; lse != NULL; prev = lse, lse = lse->next) 1.829 + if (lse->lit.col == col) break; 1.830 + xassert(lse != NULL); 1.831 + if (prev == NULL) 1.832 + set = lse->next; 1.833 + else 1.834 + prev->next = lse->next; 1.835 + dmp_free_atom(npp->pool, lse, sizeof(NPPLSE)); 1.836 + return set; 1.837 +} 1.838 + 1.839 +int npp_sat_encode_sum_ax(NPP *npp, NPPROW *row, NPPLIT y[]) 1.840 +{ NPPAIJ *aij; 1.841 + NPPLSE *set[1+NBIT_MAX], *lse; 1.842 + NPPSED sed; 1.843 + int k, n, temp; 1.844 + double sum; 1.845 + /* compute the sum (sum |a[j]|) */ 1.846 + sum = 0.0; 1.847 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.848 + sum += fabs(aij->val); 1.849 + /* determine n, the number of bits in the sum */ 1.850 + temp = (int)sum; 1.851 + if ((double)temp != sum) 1.852 + return -1; /* integer arithmetic error */ 1.853 + for (n = 0; temp > 0; n++, temp >>= 1); 1.854 + xassert(0 <= n && n <= NBIT_MAX); 1.855 + /* build initial sets J[k], 1 <= k <= n; see (5) */ 1.856 + /* set[k] is a pointer to the list of literals in J[k] */ 1.857 + for (k = 1; k <= n; k++) 1.858 + set[k] = NULL; 1.859 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.860 + { temp = (int)fabs(aij->val); 1.861 + xassert((int)temp == fabs(aij->val)); 1.862 + for (k = 1; temp > 0; k++, temp >>= 1) 1.863 + { if (temp & 1) 1.864 + { xassert(k <= n); 1.865 + lse = dmp_get_atom(npp->pool, sizeof(NPPLSE)); 1.866 + lse->lit.col = aij->col; 1.867 + lse->lit.neg = (aij->val > 0.0 ? 0 : 1); 1.868 + lse->next = set[k]; 1.869 + set[k] = lse; 1.870 + } 1.871 + } 1.872 + } 1.873 + /* main transformation loop */ 1.874 + for (k = 1; k <= n; k++) 1.875 + { /* reduce J[k] and set y[k] */ 1.876 + for (;;) 1.877 + { if (set[k] == NULL) 1.878 + { /* J[k] is empty */ 1.879 + /* set y[k] to constant false */ 1.880 + y[k].col = NULL, y[k].neg = 0; 1.881 + break; 1.882 + } 1.883 + if (set[k]->next == NULL) 1.884 + { /* J[k] contains one literal */ 1.885 + /* set y[k] to that literal */ 1.886 + y[k] = set[k]->lit; 1.887 + dmp_free_atom(npp->pool, set[k], sizeof(NPPLSE)); 1.888 + break; 1.889 + } 1.890 + if (set[k]->next->next == NULL) 1.891 + { /* J[k] contains two literals */ 1.892 + /* apply transformation (7) */ 1.893 + npp_sat_encode_sum2(npp, set[k], &sed); 1.894 + } 1.895 + else 1.896 + { /* J[k] contains at least three literals */ 1.897 + /* apply transformation (8) */ 1.898 + npp_sat_encode_sum3(npp, set[k], &sed); 1.899 + /* remove third literal from set[k] */ 1.900 + set[k] = remove_lse(npp, set[k], sed.z.col); 1.901 + } 1.902 + /* remove second literal from set[k] */ 1.903 + set[k] = remove_lse(npp, set[k], sed.y.col); 1.904 + /* remove first literal from set[k] */ 1.905 + set[k] = remove_lse(npp, set[k], sed.x.col); 1.906 + /* include new literal s to set[k] */ 1.907 + lse = dmp_get_atom(npp->pool, sizeof(NPPLSE)); 1.908 + lse->lit.col = sed.s, lse->lit.neg = 0; 1.909 + lse->next = set[k]; 1.910 + set[k] = lse; 1.911 + /* include new literal c to set[k+1] */ 1.912 + xassert(k < n); /* FIXME: can "overflow" happen? */ 1.913 + lse = dmp_get_atom(npp->pool, sizeof(NPPLSE)); 1.914 + lse->lit.col = sed.c, lse->lit.neg = 0; 1.915 + lse->next = set[k+1]; 1.916 + set[k+1] = lse; 1.917 + } 1.918 + } 1.919 + return n; 1.920 +} 1.921 + 1.922 +/*********************************************************************** 1.923 +* npp_sat_normalize_clause - normalize clause 1.924 +* 1.925 +* This routine normalizes the specified clause, which is a disjunction 1.926 +* of literals, by replacing multiple literals, which refer to the same 1.927 +* binary variable, with a single literal. 1.928 +* 1.929 +* On exit the routine returns the number of literals in the resulting 1.930 +* clause. However, if the specified clause includes both a literal and 1.931 +* its negation, the routine returns a negative value meaning that the 1.932 +* clause is equivalent to the value true. */ 1.933 + 1.934 +int npp_sat_normalize_clause(NPP *npp, int size, NPPLIT lit[]) 1.935 +{ int j, k, new_size; 1.936 + xassert(npp == npp); 1.937 + xassert(size >= 0); 1.938 + new_size = 0; 1.939 + for (k = 1; k <= size; k++) 1.940 + { for (j = 1; j <= new_size; j++) 1.941 + { if (lit[k].col == lit[j].col) 1.942 + { /* lit[k] refers to the same variable as lit[j], which 1.943 + is already included in the resulting clause */ 1.944 + if (lit[k].neg == lit[j].neg) 1.945 + { /* ignore lit[k] due to the idempotent law */ 1.946 + goto skip; 1.947 + } 1.948 + else 1.949 + { /* lit[k] is NOT lit[j]; the clause is equivalent to 1.950 + the value true */ 1.951 + return -1; 1.952 + } 1.953 + } 1.954 + } 1.955 + /* include lit[k] in the resulting clause */ 1.956 + lit[++new_size] = lit[k]; 1.957 +skip: ; 1.958 + } 1.959 + return new_size; 1.960 +} 1.961 + 1.962 +/*********************************************************************** 1.963 +* npp_sat_encode_clause - translate clause to cover inequality 1.964 +* 1.965 +* Given a clause 1.966 +* 1.967 +* OR t[j], (1) 1.968 +* j in J 1.969 +* 1.970 +* where t[j] is a literal, i.e. t[j] = x[j] or t[j] = NOT x[j], this 1.971 +* routine translates it to the following equivalent cover inequality, 1.972 +* which is added to the transformed problem: 1.973 +* 1.974 +* sum t[j] >= 1, (2) 1.975 +* j in J 1.976 +* 1.977 +* where t[j] = x[j] or t[j] = 1 - x[j]. 1.978 +* 1.979 +* If necessary, the clause should be normalized before a call to this 1.980 +* routine. */ 1.981 + 1.982 +NPPROW *npp_sat_encode_clause(NPP *npp, int size, NPPLIT lit[]) 1.983 +{ NPPROW *row; 1.984 + int k; 1.985 + xassert(size >= 1); 1.986 + row = npp_add_row(npp); 1.987 + row->lb = 1.0, row->ub = +DBL_MAX; 1.988 + for (k = 1; k <= size; k++) 1.989 + { xassert(lit[k].col != NULL); 1.990 + if (lit[k].neg == 0) 1.991 + npp_add_aij(npp, row, lit[k].col, +1.0); 1.992 + else if (lit[k].neg == 1) 1.993 + { npp_add_aij(npp, row, lit[k].col, -1.0); 1.994 + row->lb -= 1.0; 1.995 + } 1.996 + else 1.997 + xassert(lit != lit); 1.998 + } 1.999 + return row; 1.1000 +} 1.1001 + 1.1002 +/*********************************************************************** 1.1003 +* npp_sat_encode_geq - encode "not less than" constraint 1.1004 +* 1.1005 +* PURPOSE 1.1006 +* 1.1007 +* This routine translates to CNF the following constraint: 1.1008 +* 1.1009 +* n 1.1010 +* sum 2**(k-1) * y[k] >= b, (1) 1.1011 +* k=1 1.1012 +* 1.1013 +* where y[k] is either a literal (i.e. y[k] = x[k] or y[k] = 1 - x[k]) 1.1014 +* or constant false (zero), b is a given lower bound. 1.1015 +* 1.1016 +* ALGORITHM 1.1017 +* 1.1018 +* If b < 0, the constraint is redundant, so assume that b >= 0. Let 1.1019 +* 1.1020 +* n 1.1021 +* b = sum 2**(k-1) b[k], (2) 1.1022 +* k=1 1.1023 +* 1.1024 +* where b[k] is k-th binary digit of b. (Note that if b >= 2**n and 1.1025 +* therefore cannot be represented in the form (2), the constraint (1) 1.1026 +* is infeasible.) In this case the condition (1) is equivalent to the 1.1027 +* following condition: 1.1028 +* 1.1029 +* y[n] y[n-1] ... y[2] y[1] >= b[n] b[n-1] ... b[2] b[1], (3) 1.1030 +* 1.1031 +* where ">=" is understood lexicographically. 1.1032 +* 1.1033 +* Algorithmically the condition (3) can be tested as follows: 1.1034 +* 1.1035 +* for (k = n; k >= 1; k--) 1.1036 +* { if (y[k] < b[k]) 1.1037 +* y is less than b; 1.1038 +* if (y[k] > b[k]) 1.1039 +* y is greater than b; 1.1040 +* } 1.1041 +* y is equal to b; 1.1042 +* 1.1043 +* Thus, y is less than b iff there exists k, 1 <= k <= n, for which 1.1044 +* the following condition is satisfied: 1.1045 +* 1.1046 +* y[n] = b[n] AND ... AND y[k+1] = b[k+1] AND y[k] < b[k]. (4) 1.1047 +* 1.1048 +* Negating the condition (4) we have that y is not less than b iff for 1.1049 +* all k, 1 <= k <= n, the following condition is satisfied: 1.1050 +* 1.1051 +* y[n] != b[n] OR ... OR y[k+1] != b[k+1] OR y[k] >= b[k]. (5) 1.1052 +* 1.1053 +* Note that if b[k] = 0, the literal y[k] >= b[k] is always true, in 1.1054 +* which case the entire clause (5) is true and can be omitted. 1.1055 +* 1.1056 +* RETURNS 1.1057 +* 1.1058 +* Normally the routine returns zero. However, if the constraint (1) is 1.1059 +* infeasible, the routine returns non-zero. */ 1.1060 + 1.1061 +int npp_sat_encode_geq(NPP *npp, int n, NPPLIT y[], int rhs) 1.1062 +{ NPPLIT lit[1+NBIT_MAX]; 1.1063 + int j, k, size, temp, b[1+NBIT_MAX]; 1.1064 + xassert(0 <= n && n <= NBIT_MAX); 1.1065 + /* if the constraint (1) is redundant, do nothing */ 1.1066 + if (rhs < 0) 1.1067 + return 0; 1.1068 + /* determine binary digits of b according to (2) */ 1.1069 + for (k = 1, temp = rhs; k <= n; k++, temp >>= 1) 1.1070 + b[k] = temp & 1; 1.1071 + if (temp != 0) 1.1072 + { /* b >= 2**n; the constraint (1) is infeasible */ 1.1073 + return 1; 1.1074 + } 1.1075 + /* main transformation loop */ 1.1076 + for (k = 1; k <= n; k++) 1.1077 + { /* build the clause (5) for current k */ 1.1078 + size = 0; /* clause size = number of literals */ 1.1079 + /* add literal y[k] >= b[k] */ 1.1080 + if (b[k] == 0) 1.1081 + { /* b[k] = 0 -> the literal is true */ 1.1082 + goto skip; 1.1083 + } 1.1084 + else if (y[k].col == NULL) 1.1085 + { /* y[k] = 0, b[k] = 1 -> the literal is false */ 1.1086 + xassert(y[k].neg == 0); 1.1087 + } 1.1088 + else 1.1089 + { /* add literal y[k] = 1 */ 1.1090 + lit[++size] = y[k]; 1.1091 + } 1.1092 + for (j = k+1; j <= n; j++) 1.1093 + { /* add literal y[j] != b[j] */ 1.1094 + if (y[j].col == NULL) 1.1095 + { xassert(y[j].neg == 0); 1.1096 + if (b[j] == 0) 1.1097 + { /* y[j] = 0, b[j] = 0 -> the literal is false */ 1.1098 + continue; 1.1099 + } 1.1100 + else 1.1101 + { /* y[j] = 0, b[j] = 1 -> the literal is true */ 1.1102 + goto skip; 1.1103 + } 1.1104 + } 1.1105 + else 1.1106 + { lit[++size] = y[j]; 1.1107 + if (b[j] != 0) 1.1108 + lit[size].neg = 1 - lit[size].neg; 1.1109 + } 1.1110 + } 1.1111 + /* normalize the clause */ 1.1112 + size = npp_sat_normalize_clause(npp, size, lit); 1.1113 + if (size < 0) 1.1114 + { /* the clause is equivalent to the value true */ 1.1115 + goto skip; 1.1116 + } 1.1117 + if (size == 0) 1.1118 + { /* the clause is equivalent to the value false; this means 1.1119 + that the constraint (1) is infeasible */ 1.1120 + return 2; 1.1121 + } 1.1122 + /* translate the clause to corresponding cover inequality */ 1.1123 + npp_sat_encode_clause(npp, size, lit); 1.1124 +skip: ; 1.1125 + } 1.1126 + return 0; 1.1127 +} 1.1128 + 1.1129 +/*********************************************************************** 1.1130 +* npp_sat_encode_leq - encode "not greater than" constraint 1.1131 +* 1.1132 +* PURPOSE 1.1133 +* 1.1134 +* This routine translates to CNF the following constraint: 1.1135 +* 1.1136 +* n 1.1137 +* sum 2**(k-1) * y[k] <= b, (1) 1.1138 +* k=1 1.1139 +* 1.1140 +* where y[k] is either a literal (i.e. y[k] = x[k] or y[k] = 1 - x[k]) 1.1141 +* or constant false (zero), b is a given upper bound. 1.1142 +* 1.1143 +* ALGORITHM 1.1144 +* 1.1145 +* If b < 0, the constraint is infeasible, so assume that b >= 0. Let 1.1146 +* 1.1147 +* n 1.1148 +* b = sum 2**(k-1) b[k], (2) 1.1149 +* k=1 1.1150 +* 1.1151 +* where b[k] is k-th binary digit of b. (Note that if b >= 2**n and 1.1152 +* therefore cannot be represented in the form (2), the constraint (1) 1.1153 +* is redundant.) In this case the condition (1) is equivalent to the 1.1154 +* following condition: 1.1155 +* 1.1156 +* y[n] y[n-1] ... y[2] y[1] <= b[n] b[n-1] ... b[2] b[1], (3) 1.1157 +* 1.1158 +* where "<=" is understood lexicographically. 1.1159 +* 1.1160 +* Algorithmically the condition (3) can be tested as follows: 1.1161 +* 1.1162 +* for (k = n; k >= 1; k--) 1.1163 +* { if (y[k] < b[k]) 1.1164 +* y is less than b; 1.1165 +* if (y[k] > b[k]) 1.1166 +* y is greater than b; 1.1167 +* } 1.1168 +* y is equal to b; 1.1169 +* 1.1170 +* Thus, y is greater than b iff there exists k, 1 <= k <= n, for which 1.1171 +* the following condition is satisfied: 1.1172 +* 1.1173 +* y[n] = b[n] AND ... AND y[k+1] = b[k+1] AND y[k] > b[k]. (4) 1.1174 +* 1.1175 +* Negating the condition (4) we have that y is not greater than b iff 1.1176 +* for all k, 1 <= k <= n, the following condition is satisfied: 1.1177 +* 1.1178 +* y[n] != b[n] OR ... OR y[k+1] != b[k+1] OR y[k] <= b[k]. (5) 1.1179 +* 1.1180 +* Note that if b[k] = 1, the literal y[k] <= b[k] is always true, in 1.1181 +* which case the entire clause (5) is true and can be omitted. 1.1182 +* 1.1183 +* RETURNS 1.1184 +* 1.1185 +* Normally the routine returns zero. However, if the constraint (1) is 1.1186 +* infeasible, the routine returns non-zero. */ 1.1187 + 1.1188 +int npp_sat_encode_leq(NPP *npp, int n, NPPLIT y[], int rhs) 1.1189 +{ NPPLIT lit[1+NBIT_MAX]; 1.1190 + int j, k, size, temp, b[1+NBIT_MAX]; 1.1191 + xassert(0 <= n && n <= NBIT_MAX); 1.1192 + /* check if the constraint (1) is infeasible */ 1.1193 + if (rhs < 0) 1.1194 + return 1; 1.1195 + /* determine binary digits of b according to (2) */ 1.1196 + for (k = 1, temp = rhs; k <= n; k++, temp >>= 1) 1.1197 + b[k] = temp & 1; 1.1198 + if (temp != 0) 1.1199 + { /* b >= 2**n; the constraint (1) is redundant */ 1.1200 + return 0; 1.1201 + } 1.1202 + /* main transformation loop */ 1.1203 + for (k = 1; k <= n; k++) 1.1204 + { /* build the clause (5) for current k */ 1.1205 + size = 0; /* clause size = number of literals */ 1.1206 + /* add literal y[k] <= b[k] */ 1.1207 + if (b[k] == 1) 1.1208 + { /* b[k] = 1 -> the literal is true */ 1.1209 + goto skip; 1.1210 + } 1.1211 + else if (y[k].col == NULL) 1.1212 + { /* y[k] = 0, b[k] = 0 -> the literal is true */ 1.1213 + xassert(y[k].neg == 0); 1.1214 + goto skip; 1.1215 + } 1.1216 + else 1.1217 + { /* add literal y[k] = 0 */ 1.1218 + lit[++size] = y[k]; 1.1219 + lit[size].neg = 1 - lit[size].neg; 1.1220 + } 1.1221 + for (j = k+1; j <= n; j++) 1.1222 + { /* add literal y[j] != b[j] */ 1.1223 + if (y[j].col == NULL) 1.1224 + { xassert(y[j].neg == 0); 1.1225 + if (b[j] == 0) 1.1226 + { /* y[j] = 0, b[j] = 0 -> the literal is false */ 1.1227 + continue; 1.1228 + } 1.1229 + else 1.1230 + { /* y[j] = 0, b[j] = 1 -> the literal is true */ 1.1231 + goto skip; 1.1232 + } 1.1233 + } 1.1234 + else 1.1235 + { lit[++size] = y[j]; 1.1236 + if (b[j] != 0) 1.1237 + lit[size].neg = 1 - lit[size].neg; 1.1238 + } 1.1239 + } 1.1240 + /* normalize the clause */ 1.1241 + size = npp_sat_normalize_clause(npp, size, lit); 1.1242 + if (size < 0) 1.1243 + { /* the clause is equivalent to the value true */ 1.1244 + goto skip; 1.1245 + } 1.1246 + if (size == 0) 1.1247 + { /* the clause is equivalent to the value false; this means 1.1248 + that the constraint (1) is infeasible */ 1.1249 + return 2; 1.1250 + } 1.1251 + /* translate the clause to corresponding cover inequality */ 1.1252 + npp_sat_encode_clause(npp, size, lit); 1.1253 +skip: ; 1.1254 + } 1.1255 + return 0; 1.1256 +} 1.1257 + 1.1258 +/*********************************************************************** 1.1259 +* npp_sat_encode_row - encode constraint (row) of general type 1.1260 +* 1.1261 +* PURPOSE 1.1262 +* 1.1263 +* This routine translates to CNF the following constraint (row): 1.1264 +* 1.1265 +* L <= sum a[j] x[j] <= U, (1) 1.1266 +* j 1.1267 +* 1.1268 +* where all x[j] are binary variables. 1.1269 +* 1.1270 +* ALGORITHM 1.1271 +* 1.1272 +* First, the routine performs substitution x[j] = t[j] for j in J+ 1.1273 +* and x[j] = 1 - t[j] for j in J-, where J+ = { j : a[j] > 0 } and 1.1274 +* J- = { j : a[j] < 0 }. This gives: 1.1275 +* 1.1276 +* L <= sum a[j] t[j] + sum a[j] (1 - t[j]) <= U ==> 1.1277 +* j in J+ j in J- 1.1278 +* 1.1279 +* L' <= sum |a[j]| t[j] <= U', (2) 1.1280 +* j 1.1281 +* 1.1282 +* where 1.1283 +* 1.1284 +* L' = L - sum a[j], U' = U - sum a[j]. (3) 1.1285 +* j in J- j in J- 1.1286 +* 1.1287 +* (Actually only new bounds L' and U' are computed.) 1.1288 +* 1.1289 +* Then the routine translates to CNF the following equality: 1.1290 +* 1.1291 +* n 1.1292 +* sum |a[j]| t[j] = sum 2**(k-1) * y[k], (4) 1.1293 +* j k=1 1.1294 +* 1.1295 +* where y[k] is either some t[j] or a new literal or a constant zero 1.1296 +* (see the routine npp_sat_encode_sum_ax). 1.1297 +* 1.1298 +* Finally, the routine translates to CNF the following conditions: 1.1299 +* 1.1300 +* n 1.1301 +* sum 2**(k-1) * y[k] >= L' (5) 1.1302 +* k=1 1.1303 +* 1.1304 +* and 1.1305 +* 1.1306 +* n 1.1307 +* sum 2**(k-1) * y[k] <= U' (6) 1.1308 +* k=1 1.1309 +* 1.1310 +* (see the routines npp_sat_encode_geq and npp_sat_encode_leq). 1.1311 +* 1.1312 +* All resulting clauses are encoded as cover inequalities and included 1.1313 +* into the transformed problem. 1.1314 +* 1.1315 +* Note that on exit the routine removes the specified constraint (row) 1.1316 +* from the original problem. 1.1317 +* 1.1318 +* RETURNS 1.1319 +* 1.1320 +* The routine returns one of the following codes: 1.1321 +* 1.1322 +* 0 - translation was successful; 1.1323 +* 1 - constraint (1) was found infeasible; 1.1324 +* 2 - integer arithmetic error occured. */ 1.1325 + 1.1326 +int npp_sat_encode_row(NPP *npp, NPPROW *row) 1.1327 +{ NPPAIJ *aij; 1.1328 + NPPLIT y[1+NBIT_MAX]; 1.1329 + int n, rhs; 1.1330 + double lb, ub; 1.1331 + /* the row should not be free */ 1.1332 + xassert(!(row->lb == -DBL_MAX && row->ub == +DBL_MAX)); 1.1333 + /* compute new bounds L' and U' (3) */ 1.1334 + lb = row->lb; 1.1335 + ub = row->ub; 1.1336 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.1337 + { if (aij->val < 0.0) 1.1338 + { if (lb != -DBL_MAX) 1.1339 + lb -= aij->val; 1.1340 + if (ub != -DBL_MAX) 1.1341 + ub -= aij->val; 1.1342 + } 1.1343 + } 1.1344 + /* encode the equality (4) */ 1.1345 + n = npp_sat_encode_sum_ax(npp, row, y); 1.1346 + if (n < 0) 1.1347 + return 2; /* integer arithmetic error */ 1.1348 + /* encode the condition (5) */ 1.1349 + if (lb != -DBL_MAX) 1.1350 + { rhs = (int)lb; 1.1351 + if ((double)rhs != lb) 1.1352 + return 2; /* integer arithmetic error */ 1.1353 + if (npp_sat_encode_geq(npp, n, y, rhs) != 0) 1.1354 + return 1; /* original constraint is infeasible */ 1.1355 + } 1.1356 + /* encode the condition (6) */ 1.1357 + if (ub != +DBL_MAX) 1.1358 + { rhs = (int)ub; 1.1359 + if ((double)rhs != ub) 1.1360 + return 2; /* integer arithmetic error */ 1.1361 + if (npp_sat_encode_leq(npp, n, y, rhs) != 0) 1.1362 + return 1; /* original constraint is infeasible */ 1.1363 + } 1.1364 + /* remove the specified row from the problem */ 1.1365 + npp_del_row(npp, row); 1.1366 + return 0; 1.1367 +} 1.1368 + 1.1369 +/*********************************************************************** 1.1370 +* npp_sat_encode_prob - encode 0-1 feasibility problem 1.1371 +* 1.1372 +* This routine translates the specified 0-1 feasibility problem to an 1.1373 +* equivalent SAT-CNF problem. 1.1374 +* 1.1375 +* N.B. Currently this is a very crude implementation. 1.1376 +* 1.1377 +* RETURNS 1.1378 +* 1.1379 +* 0 success; 1.1380 +* 1.1381 +* GLP_ENOPFS primal/integer infeasibility detected; 1.1382 +* 1.1383 +* GLP_ERANGE integer overflow occured. */ 1.1384 + 1.1385 +int npp_sat_encode_prob(NPP *npp) 1.1386 +{ NPPROW *row, *next_row, *prev_row; 1.1387 + NPPCOL *col, *next_col; 1.1388 + int cover = 0, pack = 0, partn = 0, ret; 1.1389 + /* process and remove free rows */ 1.1390 + for (row = npp->r_head; row != NULL; row = next_row) 1.1391 + { next_row = row->next; 1.1392 + if (row->lb == -DBL_MAX && row->ub == +DBL_MAX) 1.1393 + npp_sat_free_row(npp, row); 1.1394 + } 1.1395 + /* process and remove fixed columns */ 1.1396 + for (col = npp->c_head; col != NULL; col = next_col) 1.1397 + { next_col = col->next; 1.1398 + if (col->lb == col->ub) 1.1399 + xassert(npp_sat_fixed_col(npp, col) == 0); 1.1400 + } 1.1401 + /* only binary variables should remain */ 1.1402 + for (col = npp->c_head; col != NULL; col = col->next) 1.1403 + xassert(col->is_int && col->lb == 0.0 && col->ub == 1.0); 1.1404 + /* new rows may be added to the end of the row list, so we walk 1.1405 + from the end to beginning of the list */ 1.1406 + for (row = npp->r_tail; row != NULL; row = prev_row) 1.1407 + { prev_row = row->prev; 1.1408 + /* process special cases */ 1.1409 + ret = npp_sat_is_cover_ineq(npp, row); 1.1410 + if (ret != 0) 1.1411 + { /* row is covering inequality */ 1.1412 + cover++; 1.1413 + /* since it already encodes a clause, just transform it to 1.1414 + canonical form */ 1.1415 + if (ret == 2) 1.1416 + { xassert(npp_sat_reverse_row(npp, row) == 0); 1.1417 + ret = npp_sat_is_cover_ineq(npp, row); 1.1418 + } 1.1419 + xassert(ret == 1); 1.1420 + continue; 1.1421 + } 1.1422 + ret = npp_sat_is_partn_eq(npp, row); 1.1423 + if (ret != 0) 1.1424 + { /* row is partitioning equality */ 1.1425 + NPPROW *cov; 1.1426 + NPPAIJ *aij; 1.1427 + partn++; 1.1428 + /* transform it to canonical form */ 1.1429 + if (ret == 2) 1.1430 + { xassert(npp_sat_reverse_row(npp, row) == 0); 1.1431 + ret = npp_sat_is_partn_eq(npp, row); 1.1432 + } 1.1433 + xassert(ret == 1); 1.1434 + /* and split it into covering and packing inequalities, 1.1435 + both in canonical forms */ 1.1436 + cov = npp_add_row(npp); 1.1437 + cov->lb = row->lb, cov->ub = +DBL_MAX; 1.1438 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.1439 + npp_add_aij(npp, cov, aij->col, aij->val); 1.1440 + xassert(npp_sat_is_cover_ineq(npp, cov) == 1); 1.1441 + /* the cover inequality already encodes a clause and do 1.1442 + not need any further processing */ 1.1443 + row->lb = -DBL_MAX; 1.1444 + xassert(npp_sat_is_pack_ineq(npp, row) == 1); 1.1445 + /* the packing inequality will be processed below */ 1.1446 + pack--; 1.1447 + } 1.1448 + ret = npp_sat_is_pack_ineq(npp, row); 1.1449 + if (ret != 0) 1.1450 + { /* row is packing inequality */ 1.1451 + NPPROW *rrr; 1.1452 + int nlit, desired_nlit = 4; 1.1453 + pack++; 1.1454 + /* transform it to canonical form */ 1.1455 + if (ret == 2) 1.1456 + { xassert(npp_sat_reverse_row(npp, row) == 0); 1.1457 + ret = npp_sat_is_pack_ineq(npp, row); 1.1458 + } 1.1459 + xassert(ret == 1); 1.1460 + /* process the packing inequality */ 1.1461 + for (;;) 1.1462 + { /* determine the number of literals in the remaining 1.1463 + inequality */ 1.1464 + nlit = npp_row_nnz(npp, row); 1.1465 + if (nlit <= desired_nlit) 1.1466 + break; 1.1467 + /* split the current inequality into one having not more 1.1468 + than desired_nlit literals and remaining one */ 1.1469 + rrr = npp_sat_split_pack(npp, row, desired_nlit-1); 1.1470 + /* translate the former inequality to CNF and remove it 1.1471 + from the original problem */ 1.1472 + npp_sat_encode_pack(npp, rrr); 1.1473 + } 1.1474 + /* translate the remaining inequality to CNF and remove it 1.1475 + from the original problem */ 1.1476 + npp_sat_encode_pack(npp, row); 1.1477 + continue; 1.1478 + } 1.1479 + /* translate row of general type to CNF and remove it from the 1.1480 + original problem */ 1.1481 + ret = npp_sat_encode_row(npp, row); 1.1482 + if (ret == 0) 1.1483 + ; 1.1484 + else if (ret == 1) 1.1485 + ret = GLP_ENOPFS; 1.1486 + else if (ret == 2) 1.1487 + ret = GLP_ERANGE; 1.1488 + else 1.1489 + xassert(ret != ret); 1.1490 + if (ret != 0) 1.1491 + goto done; 1.1492 + } 1.1493 + ret = 0; 1.1494 + if (cover != 0) 1.1495 + xprintf("%d covering inequalities\n", cover); 1.1496 + if (pack != 0) 1.1497 + xprintf("%d packing inequalities\n", pack); 1.1498 + if (partn != 0) 1.1499 + xprintf("%d partitioning equalities\n", partn); 1.1500 +done: return ret; 1.1501 +} 1.1502 + 1.1503 +/* eof */