001/* 002 * $Id: IntegerProgram.java 5869 2018-07-20 15:53:10Z kredel $ 003 */ 004 005package edu.jas.application; 006 007 008import java.io.IOException; 009import java.io.Reader; 010import java.io.StringReader; 011import java.util.Arrays; 012 013import org.apache.logging.log4j.Logger; 014import org.apache.logging.log4j.LogManager; 015 016import edu.jas.arith.BigInteger; 017import edu.jas.poly.ExpVector; 018import edu.jas.poly.ExpVectorLong; 019import edu.jas.poly.GenPolynomial; 020import edu.jas.poly.GenPolynomialTokenizer; 021import edu.jas.poly.PolynomialList; 022import edu.jas.poly.TermOrder; 023 024 025/** 026 * Solution of Integer Programming problems using Groebner bases. Integer 027 * Program is in standard form -> minimization of given Equation plus 028 * restrictions. See chapter 8 in Cox, Little, O'Shea "Using Algebraic 029 * Geometry", 1998. 030 * @author Maximilian Nohr 031 */ 032public class IntegerProgram implements java.io.Serializable { 033 034 035 private static final Logger logger = LogManager.getLogger(IntegerProgram.class); 036 037 038 private static boolean DEBUG = logger.isDebugEnabled(); //false; 039 040 041 private boolean negVars; 042 043 044 private boolean success; 045 046 047 /* 048 Integer Program is in standard form -> minimization of given 049 Equation + restrictions 050 */ 051 int n; // # of variables including slack variables 052 053 054 int m; // # of restrictions 055 056 057 long[] C; // List of Coefficients c_1...c_n of objective function 058 059 060 long[] B; // List of b_1...b_m, restriction right hand side 061 062 063 long[][] A; // m x n Matrix of a_{11}....a_{mn}, restriction matrix 064 065 066 long[] D; // Polynomial degrees 1...n 067 068 069 long[][] Aa; // restriction matrix a_{11}..a_{mn} after Laurent transformation 070 071 072 Ideal<BigInteger> I; //the Ideal 073 074 075 Ideal<BigInteger> GB; // the Groebner base for the ideal 076 077 078 TermOrder to; // the Term order for the GB 079 080 081 PolynomialList<BigInteger> F; // The Polynomials that generate the Ideal 082 083 084 GenPolynomial<BigInteger> S; // The Polynomial that is reduced for the Solution 085 086 087 /** 088 * Constructor. Use one instance for every new problem since solve() is not 089 * reentrant. 090 */ 091 public IntegerProgram() { 092 } 093 094 095 /** 096 * Set DEBUG flag to parameter value. 097 * @param b 098 */ 099 public void setDebug(boolean b) { 100 DEBUG = b; 101 } 102 103 104 /* 105 * Setup the Ideal corresponding to the Integer Program. 106 */ 107 @SuppressWarnings("unchecked") 108 private void createIdeal() { 109 Aa = A.clone(); 110 negVars = negVarTest(); 111 String[] w = new String[n]; 112 String[] f = new String[n]; 113 String[] z = new String[m]; 114 String[] t = new String[n]; 115 StringBuilder sb = new StringBuilder(); 116 sb.append("Int("); 117 118 if (negVars) { //A or B has negative values 119 for (int i = 1; i <= n; i++) { 120 w[i - 1] = "w" + i; 121 } 122 for (int i = 1; i <= m; i++) { 123 z[i - 1] = "z" + i; 124 } 125 for (int i = 0; i < n; i++) { 126 StringBuffer h = new StringBuffer(""); 127 long min = 0; 128 for (int j = 0; j < m; j++) { 129 if (A[j][i] < min) { 130 min = A[j][i]; 131 } 132 } 133 if (min < 0) { 134 long e = -min; 135 h.append("t^" + e + " * "); 136 for (int j = 0; j < m; j++) { 137 Aa[j][i] = A[j][i] + e; 138 h.append(z[j] + "^" + Aa[j][i] + " * "); 139 } 140 } else { 141 for (int j = 0; j < m; j++) { 142 if (A[j][i] != 0) { 143 h.append(z[j] + "^" + A[j][i] + " * "); 144 } 145 } 146 } 147 f[i] = h.substring(0, h.length() - 3).toString(); 148 } 149 setDeg(); 150 setTO(); 151 for (int i = 0; i < n; i++) { 152 t[i] = f[i] + " - " + w[i]; 153 } 154 sb.append("t"); 155 for (int i = 0; i < m; i++) { 156 sb.append(",").append(z[i]); 157 } 158 for (int i = 0; i < n; i++) { 159 sb.append(",").append(w[i]); 160 } 161 sb.append(") W "); 162 //sb.append(to.weightToString().substring(6, to.weightToString().length())); 163 sb.append(to.weightToString()); 164 sb.append(" ( ( t"); 165 for (int i = 0; i < m; i++) { 166 sb.append(" * ").append(z[i]); 167 } 168 sb.append(" - 1 )"); 169 for (int i = 0; i < n; i++) { 170 sb.append(", (").append(t[i]).append(" )"); 171 } 172 sb.append(") "); 173 174 } else { //if neither A nor B contain negative values 175 for (int i = 1; i <= n; i++) { 176 w[i - 1] = "w" + i; 177 } 178 for (int i = 1; i <= m; i++) { 179 z[i - 1] = "z" + i; 180 } 181 for (int i = 0; i < n; i++) { 182 StringBuffer h = new StringBuffer(""); 183 for (int j = 0; j < m; j++) { 184 if (A[j][i] != 0) { 185 h.append(z[j] + "^" + A[j][i] + " * "); 186 } 187 } 188 f[i] = h.substring(0, h.length() - 3).toString(); 189 } 190 setDeg(); 191 setTO(); 192 for (int i = 0; i < n; i++) { 193 t[i] = f[i] + " - " + w[i]; 194 } 195 sb.append(z[0]); 196 for (int i = 1; i < m; i++) { 197 sb.append(",").append(z[i]); 198 } 199 for (int i = 0; i < n; i++) { 200 sb.append(",").append(w[i]); 201 } 202 sb.append(") W "); 203 //sb.append(to.weightToString().substring(6, to.weightToString().length())); 204 sb.append(to.weightToString()); 205 sb.append(" ( (").append(t[0]).append(")"); 206 for (int i = 1; i < n; i++) { 207 sb.append(", (").append(t[i]).append(" )"); 208 } 209 sb.append(") "); 210 } 211 if (DEBUG) { 212 logger.debug("list=" + sb.toString()); 213 } 214 215 Reader source = new StringReader(sb.toString()); 216 GenPolynomialTokenizer parser = new GenPolynomialTokenizer(source); 217 PolynomialList<BigInteger> F = null; 218 219 try { 220 F = (PolynomialList<BigInteger>) parser.nextPolynomialSet(); 221 } catch (ClassCastException e) { 222 e.printStackTrace(); 223 return; 224 } catch (IOException e) { 225 e.printStackTrace(); 226 return; 227 } 228 if (DEBUG) { 229 logger.debug("F=" + F); 230 } 231 I = new Ideal<BigInteger>(F); 232 return; 233 } 234 235 236 /** 237 * @return true if the last calculation had a solution, else false 238 */ 239 public boolean getSuccess() { 240 return success; 241 } 242 243 244 /** 245 * Solve Integer Program. 246 * @param A matrix of restrictions 247 * @param B restrictions right hand side 248 * @param C objective function 249 * @return solution s, such that s*C -> minimal and A*s = B, or s = null 250 * if no solution exists 251 */ 252 public long[] solve(long[][] A, long[] B, long[] C) { 253 this.A = Arrays.copyOf(A, A.length); 254 this.B = Arrays.copyOf(B, B.length); 255 this.C = Arrays.copyOf(C, C.length); 256 this.n = A[0].length; 257 this.m = A.length; 258 D = new long[n]; 259 260 createIdeal(); 261 GB = I.GB(); 262 return solve(B); 263 } 264 265 266 /** 267 * Solve Integer Program for new right hand side. Uses the GB (matrix A and 268 * C) from the last calculation. 269 * @param B restrictions right hand side 270 * @return solution s, such that s*C -> minimal and A*s = B, or s = null 271 * if no solution exists 272 */ 273 public long[] solve(long[] B) { 274 long[] returnMe = new long[n]; 275 if (B.length != m) { 276 System.out.println("ERROR: Dimensions don't match: " + B.length + " != " + m); 277 return returnMe; 278 } 279 long[] l; 280 this.B = Arrays.copyOf(B, B.length); 281 if (DEBUG) { 282 logger.debug("GB=" + GB); 283 } 284 if (negVars) { 285 l = new long[m + n + 1]; 286 long min = findMin(B); 287 if (min < 0) { 288 long r = -min; 289 l[m + n] = r; 290 for (int i = 0; i < m; i++) { 291 l[m + n - 1 - i] = B[i] + r; 292 } 293 } else { 294 for (int i = 0; i < m; i++) { 295 l[m + n - 1 - i] = B[i]; 296 } 297 } 298 } else { 299 l = new long[m + n]; 300 for (int i = 0; i < m; i++) { 301 l[m + n - 1 - i] = B[i]; 302 } 303 } 304 ExpVector e = new ExpVectorLong(l); 305 S = new GenPolynomial<BigInteger>(I.getRing(), e); 306 S = GB.normalform(S); 307 308 ExpVector E = S.exponentIterator().next(); 309 for (int i = 0; i < n; i++) { 310 returnMe[n - 1 - i] = E.getVal(i); 311 } 312 success = true; 313 for (int i = n; i < n + m; i++) { 314 if (E.getVal(i) != 0) { 315 success = false; 316 break; 317 } 318 } 319 if (success) { 320 if (DEBUG) { 321 logger.debug("The solution is: " + Arrays.toString(returnMe)); 322 } 323 } else { 324 logger.warn("The Problem does not have a feasible solution."); 325 returnMe = null; 326 } 327 return returnMe; 328 } 329 330 331 /* 332 * Set the degree. 333 */ 334 private void setDeg() { 335 for (int j = 0; j < n; j++) { 336 for (int i = 0; i < m; i++) { 337 D[j] += Aa[i][j]; 338 } 339 } 340 } 341 342 343 /* 344 * Set the term order. 345 */ 346 private void setTO() { 347 int h; 348 if (negVars) {//if A and/or B contains negative values 349 h = m + n + 1; 350 } else { 351 h = m + n; 352 } 353 long[] u1 = new long[h]; 354 long[] u2 = new long[h]; 355 for (int i = 0; i < h - n; i++) { //m+1 because t needs another 1 356 u1[h - 1 - i] = 1; 357 } 358 long[] h1 = new long[h]; // help vectors to construct u2 out of 359 long[] h2 = new long[h]; 360 for (int i = h - n; i < h; i++) { 361 h1[i] = C[i - (h - n)]; 362 h2[i] = D[i - (h - n)]; 363 } 364 long min = h1[0]; 365 for (int i = 0; i < h; i++) { 366 u2[h - 1 - i] = h1[i] + h2[i]; 367 if (u2[h - 1 - i] < min) { 368 min = u2[h - 1 - i]; 369 } 370 } 371 while (min < 0) { 372 min = u2[0]; 373 for (int i = 0; i < h; i++) { 374 u2[h - 1 - i] += h2[i]; 375 if (u2[h - 1 - i] < min) { 376 min = u2[h - 1 - i]; 377 } 378 } 379 } 380 long[][] wv = { u1, u2 }; 381 to = new TermOrder(wv); 382 } 383 384 385 /* 386 * @see java.lang.Object#toString() 387 */ 388 @Override 389 public String toString() { 390 StringBuilder sb = new StringBuilder(); 391 sb.append("Function to minimize:\n"); 392 393 char c = 'A'; // variables are named A, B, C, ... 394 395 boolean plus = false; 396 for (int i = 0; i < n; i++) { 397 if (C[i] != 0) { 398 if (C[i] < 0) { 399 sb.append("(").append(C[i]).append(")"); 400 sb.append("*"); 401 } else if (C[i] != 1) { 402 sb.append(C[i]); 403 sb.append("*"); 404 } 405 sb.append(c); 406 sb.append(" + "); 407 plus = true; 408 } 409 c++; 410 } 411 if (plus) { 412 sb.delete(sb.lastIndexOf("+"), sb.length()); 413 } 414 sb.append("\nunder the Restrictions:\n"); 415 for (int i = 0; i < m; i++) { 416 c = 'A'; 417 //System.out.println("A["+i+"] = " + Arrays.toString(A[i])); 418 plus = false; 419 for (int j = 0; j < n; j++) { 420 if (A[i][j] != 0) { 421 if (A[i][j] < 0) { 422 sb.append("(").append(A[i][j]).append(")"); 423 sb.append("*"); 424 } else if (A[i][j] != 1) { 425 sb.append(A[i][j]); 426 sb.append("*"); 427 } 428 sb.append(c); 429 sb.append(" + "); 430 plus = true; 431 } 432 c++; 433 } 434 if (plus) { 435 sb.delete(sb.lastIndexOf("+"), sb.length()); 436 } else { 437 sb.append("0 "); 438 } 439 sb.append("= ").append(B[i]).append("\n"); 440 } 441 return sb.toString(); 442 } 443 444 445 /* 446 * Test for negative variables. 447 * @return true if negative variables appear 448 */ 449 private boolean negVarTest() { 450 for (int i = 0; i < m; i++) { 451 if (B[i] < 0) { 452 return true; 453 } 454 for (int j = 0; j < n; j++) { 455 if (A[i][j] < 0) { 456 return true; 457 } 458 459 } 460 } 461 return false; 462 } 463 464 465 /* 466 * Find minimal element. 467 * @param B vector of at least one element, B.length >= 1 468 * @return minimal element of B 469 */ 470 private long findMin(long[] B) { 471 long min = B[0]; 472 for (int i = 1; i < B.length; i++) { 473 if (B[i] < min) { 474 min = B[i]; 475 } 476 } 477 return min; 478 } 479 480}