Source code

001/*
002 * $Id$
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 -&gt; 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 -&gt; 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);
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 -&gt; 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 -&gt; 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)); // or () -> toString(ll)
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 &gt;= 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}