diff options
| author |  tabarbe <unknown> | 2010-08-26 23:37:01 +0000 |
|---|---|---|
| committer | tabarbe <unknown> | 2010-08-26 23:37:01 +0000 |
| commit | 47171ab9f9d31dab0d5e0a4c3c95c763452e9295 (patch) | |
| tree | 402d453ee1c63dff1a04d03eabfc2bef32eed4ed /Source/AIFramework/Polyhedra/SimplexTableau.cs | |
| parent | 8b0392fe672ce820ba07af673fe9177babdee00b (diff) | |
Boogie: Renaming the AIFramework sources in preparation for committal of my port of the project
Diffstat (limited to 'Source/AIFramework/Polyhedra/SimplexTableau.cs')
| -rw-r--r-- | Source/AIFramework/Polyhedra/SimplexTableau.cs | 717 |
1 files changed, 717 insertions, 0 deletions
diff --git a/Source/AIFramework/Polyhedra/SimplexTableau.cs b/Source/AIFramework/Polyhedra/SimplexTableau.cs new file mode 100644 index 00000000..b6f4095c --- /dev/null +++ b/Source/AIFramework/Polyhedra/SimplexTableau.cs @@ -0,0 +1,717 @@ +//-----------------------------------------------------------------------------
+//
+// Copyright (C) Microsoft Corporation. All Rights Reserved.
+//
+//-----------------------------------------------------------------------------
+namespace Microsoft.AbstractInterpretationFramework
+{
+ using System.Collections;
+ using System;
+ using Microsoft.Contracts;
+ using Microsoft.Basetypes;
+ using IMutableSet = Microsoft.Boogie.Set;
+ using HashSet = Microsoft.Boogie.Set;
+
+
+ /// <summary>
+ /// Used by LinearConstraintSystem.GenerateFrameFromConstraints.
+ /// </summary>
+ public class SimplexTableau
+ {
+ readonly int rows;
+ readonly int columns;
+ readonly Rational[,]! m;
+
+ readonly int numInitialVars;
+ readonly int numSlackVars;
+ readonly int rhsColumn;
+
+ readonly ArrayList /*IVariable!*/! dims;
+ readonly int[]! basisColumns;
+ readonly int[]! inBasis;
+ bool constructionDone = false;
+
+ void CheckInvariant()
+ {
+ assert(rows == m.GetLength(0));
+ assert(1 <= columns && columns == m.GetLength(1));
+ assert(0 <= numInitialVars);
+ assert(0 <= numSlackVars && numSlackVars <= rows);
+ assert(numInitialVars + numSlackVars + 1 == columns);
+ assert(rhsColumn == columns - 1);
+ assert(dims.Count == numInitialVars);
+ assert(basisColumns.Length == rows);
+ assert(inBasis.Length == numInitialVars + numSlackVars);
+
+ bool[] b = new bool[numInitialVars + numSlackVars];
+ int numColumnsInBasis = 0;
+ int numUninitializedRowInfo = 0;
+ for (int i = 0; i < rows; i++)
+ {
+ int c = basisColumns[i];
+ if (c == rhsColumn)
+ {
+ // all coefficients in this row are 0 (but the right-hand side may be non-0)
+ for (int j = 0; j < rhsColumn; j++)
+ {
+ assert m[i,j].IsZero;
+ }
+ numColumnsInBasis++;
+ }
+ else if (c == -1)
+ {
+ assert(!constructionDone);
+ numUninitializedRowInfo++;
+ }
+ else
+ {
+ // basis column is a column
+ assert(0 <= c && c < numInitialVars + numSlackVars);
+ // basis column is unique
+ assert(!b[c]);
+ b[c] = true;
+ // column is marked as being in basis
+ assert(inBasis[c] == i);
+ // basis column really is a basis column
+ for (int j = 0; j < rows; j++)
+ {
+ if (j == i)
+ {
+ assert m[j,c].HasValue(1);// == (Rational)new Rational(1));
+ }
+ else
+ {
+ assert m[j,c].IsZero;
+ }
+ }
+ }
+ }
+ // no other columns are marked as being in basis
+ foreach (int i in inBasis)
+ {
+ if (0 <= i)
+ {
+ assert(i < rows);
+ numColumnsInBasis++;
+ }
+ else
+ {
+ assert(i == -1);
+ }
+ }
+ assert(rows - numUninitializedRowInfo <= numColumnsInBasis && numColumnsInBasis <= rows);
+ assert(!constructionDone || numUninitializedRowInfo == 0);
+ }
+
+ /// <summary>
+ /// Constructs a matrix that represents the constraints "constraints", adding slack
+ /// variables for the inequalities among "constraints". Puts the matrix in canonical
+ /// form.
+ /// </summary>
+ /// <param name="constraints"></param>
+ [NotDelayed]
+ public SimplexTableau(ArrayList /*LinearConstraint*/! constraints)
+ {
+#if DEBUG_PRINT
+ Console.WriteLine("DEBUG: SimplexTableau constructor called with:");
+ foreach (LinearConstraint lc in constraints)
+ {
+ Console.WriteLine(" {0}", lc);
+ }
+#endif
+ // Note: This implementation is not particularly efficient, but it'll do for now.
+
+ ArrayList dims = this.dims = new ArrayList /*IVariable!*/ ();
+ int slacks = 0;
+ foreach (LinearConstraint! cc in constraints)
+ {
+ foreach (IVariable! dim in cc.coefficients.Keys)
+ {
+ if (!dims.Contains(dim))
+ {
+ dims.Add(dim);
+ }
+ }
+ if (cc.Relation == LinearConstraint.ConstraintRelation.LE)
+ {
+ slacks++;
+ }
+ }
+
+ int numInitialVars = this.numInitialVars = dims.Count;
+ int numSlackVars = this.numSlackVars = slacks;
+ int rows = this.rows = constraints.Count;
+ int columns = this.columns = numInitialVars + numSlackVars + 1;
+ this.m = new Rational[rows, columns];
+ this.rhsColumn = columns-1;
+ this.basisColumns = new int[rows];
+ this.inBasis = new int[columns-1];
+
+ base();
+
+ for (int i = 0; i < inBasis.Length; i++)
+ {
+ inBasis[i] = -1;
+ }
+
+ // Fill in the matrix
+ int r = 0;
+ int iSlack = 0;
+ foreach (LinearConstraint! cc in constraints)
+ {
+ for (int i = 0; i < dims.Count; i++)
+ {
+ m[r,i] = cc[(IVariable!)dims[i]];
+ }
+ if (cc.Relation == LinearConstraint.ConstraintRelation.LE)
+ {
+ m[r, numInitialVars + iSlack] = Rational.ONE;
+ basisColumns[r] = numInitialVars + iSlack;
+ inBasis[numInitialVars + iSlack] = r;
+ iSlack++;
+ }
+ else
+ {
+ basisColumns[r] = -1; // special value to communicate to Pivot that basis column i hasn't been set up yet
+ }
+ m[r,rhsColumn] = cc.rhs;
+ r++;
+ }
+ assert(r == constraints.Count);
+ assert(iSlack == numSlackVars);
+#if DEBUG_PRINT
+ Console.WriteLine("DEBUG: Intermediate tableau state in SimplexTableau constructor:");
+ Dump();
+#endif
+
+ // Go through the rows with uninitialized basis columns. These correspond to equality constraints.
+ // For each one, find an initial variable (non-slack variable) whose column we can make the basis
+ // column of the row.
+ for (int i = 0; i < rows; i++)
+ {
+ if (basisColumns[i] != -1)
+ {
+ continue;
+ }
+ // Find a non-0 column in row i that we can make a basis column. Since rows corresponding
+ // to equality constraints don't have slack variables and since the pivot operations performed
+ // by iterations of this loop don't introduce any non-0 coefficients in the slack-variable
+ // columns of these rows, we only need to look through the columns corresponding to initial
+ // variables.
+ for (int j = 0; j < numInitialVars; j++)
+ {
+ if (m[i,j].IsNonZero)
+ {
+#if DEBUG_PRINT
+ Console.WriteLine("-- About to Pivot({0},{1})", i, j);
+#endif
+ assert(inBasis[j] == -1);
+ Pivot(i,j);
+#if DEBUG_PRINT
+ Console.WriteLine("Tableau after Pivot:");
+ Dump();
+ #endif
+ goto SET_UP_NEXT_INBASIS_COLUMN;
+ }
+ }
+ // Check the assertion in the comment above, that is, that columns corresponding to slack variables
+ // are 0 in this row.
+ for (int j = numInitialVars; j < rhsColumn; j++)
+ {
+ assert m[i,j].IsZero;
+ }
+ // There is no column in this row that we can put into basis.
+ basisColumns[i] = rhsColumn;
+ SET_UP_NEXT_INBASIS_COLUMN: {}
+ }
+
+ constructionDone = true;
+ CheckInvariant();
+ }
+
+ public IMutableSet! /*IVariable!*/ GetDimensions()
+ {
+ HashSet /*IVariable!*/ z = new HashSet /*IVariable!*/ ();
+ foreach (IVariable! dim in dims)
+ {
+ z.Add(dim);
+ }
+ return z;
+ }
+
+ public Rational this [int r, int c]
+ {
+ get
+ {
+ return m[r,c];
+ }
+ set
+ {
+ m[r,c] = value;
+ }
+ }
+
+ /// <summary>
+ /// Applies the Pivot Operation on row "r" and column "c".
+ ///
+ /// This method can be called when !constructionDone, that is, at a time when not all basis
+ /// columns have been set up (indicated by -1 in basisColumns). This method helps set up
+ /// those basis columns.
+ ///
+ /// The return value is an undo record that can be used with UnPivot.
+ /// </summary>
+ /// <param name="r"></param>
+ /// <param name="c"></param>
+ public Rational[]! Pivot(int r, int c)
+ {
+ assert(0 <= r && r < rows);
+ assert(0 <= c && c < columns-1);
+ assert(m[r,c].IsNonZero);
+ assert(inBasis[c] == -1); // follows from invariant and m[r,c] != 0
+ assert(basisColumns[r] != rhsColumn); // follows from invariant and m[r,c] != 0
+
+ Rational[] undo = new Rational[rows+1];
+ for (int i = 0; i < rows; i++)
+ {
+ undo[i] = m[i,c];
+ }
+
+ // scale the pivot row
+ Rational q = m[r,c];
+ if (q != Rational.ONE)
+ {
+ for (int j = 0; j < columns; j++)
+ {
+ m[r,j] /= q;
+ }
+ }
+
+ // subtract a multiple of the pivot row from all other rows
+ for (int i = 0; i < rows; i++)
+ {
+ if (i != r)
+ {
+ q = m[i,c];
+ if (q.IsNonZero)
+ {
+ for (int j = 0; j < columns; j++)
+ {
+ m[i,j] -= q * m[r,j];
+ }
+ }
+ }
+ }
+
+ // update basis information
+ int prevCol = basisColumns[r];
+ undo[rows] = Rational.FromInt(prevCol);
+ basisColumns[r] = c;
+ if (prevCol != -1)
+ {
+ inBasis[prevCol] = -1;
+ }
+ inBasis[c] = r;
+
+ return undo;
+ }
+
+ /// <summary>
+ /// If the last operation applied to the tableau was:
+ /// undo = Pivot(i,j);
+ /// then UnPivot(i, j, undo) undoes the pivot operation.
+ /// Note: This operation is not supported for any call to Pivot before constructionDone
+ /// is set to true.
+ /// </summary>
+ /// <param name="r"></param>
+ /// <param name="c"></param>
+ /// <param name="undo"></param>
+ void UnPivot(int r, int c, Rational[]! undo)
+ {
+ assert(0 <= r && r < rows);
+ assert(0 <= c && c < columns-1);
+ assert(m[r,c].HasValue(1));
+ assert(undo.Length == rows+1);
+
+ // add a multiple of the pivot row to all other rows
+ for (int i = 0; i < rows; i++)
+ {
+ if (i != r)
+ {
+ Rational q = undo[i];
+ if (q.IsNonZero)
+ {
+ for (int j = 0; j < columns; j++)
+ {
+ m[i,j] += q * m[r,j];
+ }
+ }
+ }
+ }
+
+ // scale the pivot row
+ Rational p = undo[r];
+ for (int j = 0; j < columns; j++)
+ {
+ m[r,j] *= p;
+ }
+
+ // update basis information
+ int prevCol = undo[rows].AsInteger;
+ assert(prevCol != -1);
+ basisColumns[r] = prevCol;
+ inBasis[c] = -1;
+ inBasis[prevCol] = r;
+ }
+
+ /// <summary>
+ /// Returns true iff the current basis of the system of constraints modeled by the simplex tableau
+ /// is feasible. May have a side effect of performing a number of pivot operations on the tableau,
+ /// but any such pivot operation will be in the columns of slack variables (that is, this routine
+ /// does not change the set of initial-variable columns in basis).
+ ///
+ /// CAVEAT: I have no particular reason to believe that the algorithm used here will terminate. --KRML
+ /// </summary>
+ /// <returns></returns>
+ public bool IsFeasibleBasis
+ {
+ get
+ {
+ // while there is a slack variable in basis whose row has a negative right-hand side
+ while (true)
+ {
+ bool feasibleBasis = true;
+ for (int c = numInitialVars; c < rhsColumn; c++)
+ {
+ int k = inBasis[c];
+ if (0 <= k && k < rhsColumn && m[k,rhsColumn].IsNegative)
+ {
+ assert(m[k,c].HasValue(1)); // c is in basis
+ // Try to pivot on a different slack variable in this row
+ for (int i = numInitialVars; i < rhsColumn; i++)
+ {
+ if (m[k,i].IsNegative)
+ {
+ assert(c != i); // c is in basis, so m[k,c]==1, which is not negative
+ Pivot(k, i);
+#if DEBUG_PRINT
+ Console.WriteLine("Tableau after Pivot operation on ({0},{1}) in IsFeasibleBasis:", k, i);
+ Dump();
+#endif
+ assert(inBasis[c] == -1);
+ assert(inBasis[i] == k);
+ assert(m[k,rhsColumn].IsNonNegative);
+ goto START_ANEW;
+ }
+ }
+ feasibleBasis = false;
+ }
+ }
+ return feasibleBasis;
+ START_ANEW: ;
+ }
+ return false; // make compiler shut up
+ }
+ }
+
+ /// <summary>
+ /// Whether or not all initial variables (the non-slack variables) are in basis)
+ /// </summary>
+ public bool AllInitialVarsInBasis
+ {
+ get
+ {
+ for (int i = 0; i < numInitialVars; i++)
+ {
+ if (inBasis[i] == -1)
+ {
+ return false;
+ }
+ }
+ return true;
+ }
+ }
+
+ /// <summary>
+ /// Adds as many initial variables as possible to the basis.
+ /// </summary>
+ /// <returns></returns>
+ public void AddInitialVarsToBasis()
+ {
+ // while there exists an initial variable not in the basis and not satisfying
+ // condition 3.4.2.2 in Cousot and Halbwachs, perform a pivot operation
+ while (true)
+ {
+ for (int i = 0; i < numInitialVars; i++)
+ {
+ if (inBasis[i] == -1)
+ {
+ // initial variable i is not in the basis
+ for (int j = 0; j < rows; j++)
+ {
+ if (m[j,i].IsNonZero)
+ {
+ int k = basisColumns[j];
+ if (numInitialVars <= k && k < rhsColumn)
+ {
+ // slack variable k is in basis for row j
+ Pivot(j, i);
+ assert(inBasis[k] == -1);
+ assert(inBasis[i] == j && basisColumns[j] == i);
+ goto START_ANEW;
+ }
+ }
+ }
+ }
+ }
+ // No more initial variables can be moved into basis.
+ return;
+ START_ANEW: {}
+ }
+ }
+
+ /// <summary>
+ /// Adds to "lines" the lines implied by initial-variable columns not in basis
+ /// (see section 3.4.2 of Cousot and Halbwachs), and adds to "constraints" the
+ /// constraints to exclude those lines (see step 4.2 of section 3.4.3 of
+ /// Cousot and Halbwachs).
+ /// </summary>
+ /// <param name="lines"></param>
+ /// <param name="constraints"></param>
+ public void ProduceLines(ArrayList /*FrameElement*/! lines, ArrayList /*LinearConstraint*/! constraints)
+ {
+ // for every initial variable not in basis
+ for (int i0 = 0; i0 < numInitialVars; i0++)
+ {
+ if (inBasis[i0] == -1)
+ {
+ FrameElement fe = new FrameElement();
+ LinearConstraint lc = new LinearConstraint(LinearConstraint.ConstraintRelation.EQ);
+ for (int i = 0; i < numInitialVars; i++)
+ {
+ if (i == i0)
+ {
+ fe.AddCoordinate((IVariable!)dims[i], Rational.ONE);
+ lc.SetCoefficient((IVariable!)dims[i], Rational.ONE);
+ }
+ else if (inBasis[i] != -1)
+ {
+ // i is a basis column
+ assert(m[inBasis[i],i].HasValue(1));
+ Rational val = -m[inBasis[i],i0];
+ fe.AddCoordinate((IVariable!)dims[i], val);
+ lc.SetCoefficient((IVariable!)dims[i], val);
+ }
+ }
+ lines.Add(fe);
+ constraints.Add(lc);
+ }
+ }
+ }
+
+ /// <summary>
+ /// From a feasible point where all initial variables are in the basis, traverses
+ /// all feasible bases containing all initial variables. For each such basis, adds
+ /// the vertices to "vertices" and adds to "rays" the extreme rays. See step 4.2
+ /// in section 3.4.3 of Cousot and Halbwachs.
+ /// A more efficient algorithm is found in the paper "An algorithm for
+ /// determining all extreme points of a convex polytope" by N. E. Dyer and L. G. Proll,
+ /// Mathematical Programming, 12, 1977.
+ /// Assumes that the tableau is in a state where all initial variables are in the basis.
+ /// This method has no net effect on the tableau.
+ /// Note: Duplicate vertices and rays may be added.
+ /// </summary>
+ /// <param name="vertices"></param>
+ /// <param name="rays"></param>
+ public void TraverseVertices(ArrayList! /*FrameElement*/ vertices, ArrayList! /*FrameElement*/ rays)
+ {
+ ArrayList /*bool[]*/ basesSeenSoFar = new ArrayList /*bool[]*/ ();
+ TraverseBases(basesSeenSoFar, vertices, rays);
+ }
+
+ /// <summary>
+ /// Worker method of TraverseVertices.
+ /// This method has no net effect on the tableau.
+ /// </summary>
+ /// <param name="basesSeenSoFar"></param>
+ /// <param name="vertices"></param>
+ /// <param name="rays"></param>
+ void TraverseBases(ArrayList /*bool[]*/! basesSeenSoFar, ArrayList /*FrameElement*/! vertices, ArrayList /*FrameElement*/! rays)
+ {
+ CheckInvariant();
+
+ bool[] thisBasis = new bool[numSlackVars];
+ for (int i = numInitialVars; i < rhsColumn; i++)
+ {
+ if (inBasis[i] != -1)
+ {
+ thisBasis[i-numInitialVars] = true;
+ }
+ }
+ foreach (bool[]! basis in basesSeenSoFar)
+ {
+ assert(basis.Length == numSlackVars);
+ for (int i = 0; i < numSlackVars; i++)
+ {
+ if (basis[i] != thisBasis[i])
+ {
+ goto COMPARE_WITH_NEXT_BASIS;
+ }
+ }
+ // thisBasis and basis are the same--that is, basisColumns has been visited before--so
+ // we don't traverse anything from here
+ return;
+ COMPARE_WITH_NEXT_BASIS: {}
+ }
+ // basisColumns has not been seen before; record thisBasis and continue with the traversal here
+ basesSeenSoFar.Add(thisBasis);
+
+#if DEBUG_PRINT
+ Console.Write("TraverseBases, new basis: ");
+ foreach (bool t in thisBasis) {
+ Console.Write("{0}", t ? "*" : ".");
+ }
+ Console.WriteLine();
+ Dump();
+#endif
+ // Add vertex
+ FrameElement v = new FrameElement();
+ for (int i = 0; i < rows; i++)
+ {
+ int j = basisColumns[i];
+ if (j < numInitialVars)
+ {
+ v.AddCoordinate((IVariable!)dims[j], m[i,rhsColumn]);
+ }
+ }
+#if DEBUG_PRINT
+ Console.WriteLine(" Adding vertex: {0}", v);
+#endif
+ vertices.Add(v);
+
+ // Add rays. Traverse all columns corresponding to slack variables that
+ // are not in basis (see second bullet of section 3.4.2 of Cousot and Halbwachs).
+ for (int i0 = numInitialVars; i0 < rhsColumn; i0++)
+ {
+ if (inBasis[i0] != -1)
+ {
+ // skip those slack-variable columns that are in basis
+ continue;
+ }
+ // check if slack-variable, non-basis column i corresponds to an extreme ray
+ for (int row = 0; row < rows; row++)
+ {
+ if (m[row,i0].IsPositive)
+ {
+ for (int k = numInitialVars; k < rhsColumn; k++)
+ {
+ if (inBasis[k] != -1 && m[row,k].IsNonZero)
+ {
+ // does not correspond to an extreme ray
+ goto CHECK_NEXT_SLACK_VAR;
+ }
+ }
+ }
+ }
+ // corresponds to an extreme ray
+ FrameElement ray = new FrameElement();
+ for (int i = 0; i < numInitialVars; i++)
+ {
+ int j0 = inBasis[i];
+ Rational val = -m[j0,i0];
+ ray.AddCoordinate((IVariable!)dims[i], val);
+ }
+#if DEBUG_PRINT
+ Console.WriteLine(" Adding ray: {0}", ray);
+#endif
+ rays.Add(ray);
+ CHECK_NEXT_SLACK_VAR: {}
+ }
+
+ // Continue traversal
+ for (int i = numInitialVars; i < rhsColumn; i++)
+ {
+ int j = inBasis[i];
+ if (j != -1)
+ {
+ // try moving i out of basis and some other slack-variable column into basis
+ for (int k = numInitialVars; k < rhsColumn; k++)
+ {
+ if (inBasis[k] == -1 && m[j,k].IsPositive)
+ {
+ Rational[] undo = Pivot(j, k);
+ // check if the new basis is feasible
+ for (int p = 0; p < rows; p++)
+ {
+ int c = basisColumns[p];
+ if (numInitialVars <= c && c < rhsColumn && m[p,rhsColumn].IsNegative)
+ {
+ // not feasible
+ goto AFTER_TRAVERSE;
+ }
+ }
+ TraverseBases(basesSeenSoFar, vertices, rays);
+ AFTER_TRAVERSE:
+ UnPivot(j, k, undo);
+ }
+ }
+ }
+ }
+ }
+
+ public void Dump()
+ {
+ // names
+ Console.Write(" ");
+ for (int i = 0; i < numInitialVars; i++)
+ {
+ Console.Write(" {0,4} ", dims[i]);
+ }
+ Console.WriteLine();
+ // numbers
+ Console.Write(" ");
+ for (int i = 0; i < columns; i++)
+ {
+ if (i == numInitialVars || i == rhsColumn)
+ {
+ Console.Write("|");
+ }
+ Console.Write(" {0,4}", i);
+ if (i < rhsColumn && inBasis[i] != -1)
+ {
+ Console.Write("* ");
+ assert(basisColumns[inBasis[i]] == i);
+ }
+ else
+ {
+ Console.Write(" ");
+ }
+ }
+ Console.WriteLine();
+ // line
+ Console.Write(" ");
+ for (int i = 0; i < columns; i++)
+ {
+ if (i == numInitialVars || i == rhsColumn)
+ {
+ Console.Write("+");
+ }
+ Console.Write("---------");
+ }
+ Console.WriteLine();
+
+ for (int j = 0; j < rows; j++)
+ {
+ Console.Write("{0,4}: ", basisColumns[j]);
+ for (int i = 0; i < columns; i++)
+ {
+ if (i == numInitialVars || i == rhsColumn)
+ {
+ Console.Write("|");
+ }
+ Console.Write(" {0,4:n1} ", m[j,i]);
+ }
+ Console.WriteLine();
+ }
+ }
+ }
+}
|