[sympycore] r1228 committed - Implemented crisscross algorithm for solving LP problems.
codesite...@google.com
Author: pearu.peterson
Date: Sun Jan 15 14:44:56 2012
Log: Implemented crisscross algorithm for solving LP problems.
http://code.google.com/p/sympycore/source/detail?r=1228
Added:
/trunk/sympycore/matrices/linalg_lp.py
Modified:
/trunk/sympycore/__init__.py
/trunk/sympycore/matrices/algebra.py
/trunk/sympycore/matrices/linalg.py
/trunk/sympycore/matrices/polyhedra.py
/trunk/sympycore/matrices/tests/test_polyhedra.py
=======================================
--- /dev/null
+++ /trunk/sympycore/matrices/linalg_lp.py Sun Jan 15 14:44:56 2012
@@ -0,0 +1,155 @@
+# Author: Pearu Peterson
+# Created: January 2012
+
+import copy
+from sympycore.arithmetic.numbers import div
+from .linalg import get_rc_maps
+from .algebra import Matrix
+
+class LPError(Exception):
+ """Generic Python-exception-derived object raised by sympycore.linalg
LP related functions.
+ """
+ pass
+
+def MATRIX_DICT_LP_solve(self, method='crisscross', overwrite=False):
+ """ Solve LP problem ``max(c^T x: A*x <= b & x>=0)``.
+
+ Given matrix must contain the problem data in dictionary form::
+
+ [ 0 c^T ]
+ D = [ ]
+ [ b -A ]
+
+ Parameters
+ ----------
+ method : {'crisscross'}
+ Specify the method to be used for solving the LP problem.
+
+ overwrite : bool
+ When True then discard the content of given matrix. After
+ computation the dictionary matrix will be in either feasible or
+ inconsistent form.
+
+ Returns
+ -------
+ xopt : Matrix
+ Basic solution of the LP problem when the LP problem is feasible.
+ vopt : number
+ Optimal value of the LP problem: ``c^T xopt``
+
+
+ When the LP problem is inconsistent, LPError exception will be
+ returned.
+ """
+ head, data = self.pair
+ m, n = head.shape
+ if not (overwrite and self.is_writable):
+ data = dict(data)
+ if method=='crisscross':
+ if head.is_transpose or head.is_diagonal:
+ raise NotImplementedError('LP_solve_crisscross_MATRIX_T,
LP_solve_crisscross_MATRIX_D')
+ return LP_solve_crisscross_MATRIX(m, n, data)
+ else:
+ raise NotImplementedError(`method`)
+
+def LP_solve_crisscross_MATRIX(m, n, data):
+ """ Solve LP problem using crisscross method. See
+ MATRIX_DICT_LP_solve for definition and user interface.
+
+ The implementation of this function is based on the crisscross
+ algorithm given in:
+
+
http://www.ifor.math.ethz.ch/teaching/Courses/Fall_2011/intro_fall_11/Script2
+
+ The algorithm has few modifications to achieve more optimal Python
+ implementation. For instance, the code uses native indexing of the
+ dictionary matrix instead of label indexing. Also, finding the
+ minimum of labels is replaced with finding the maximum of indexes.
+
+ The current implementation is compact, that is, all pivot
+ operations are applied in-situ on an input matrix taking full
+ advantage of possible sparsity of the input matrix.
+ """
+ N = range(1,n) # variables are indices
+ B = range(-1,-m,-1) # slack variables are negative indices
+
+ while True:
+ Drows, Dcols = get_rc_maps(data)
+
+ def rowmax(j, Dcols, B):
+ k, index = None, None
+ for i in Dcols.get(j,[]):
+ if i and data[i, j] < 0:
+ label = B[i-1]
+ if label > k:
+ k = label
+ index = i
+ return k, index
+ def colmax(i, Drows, N):
+ k, index = None, None
+ for j in Drows.get(i,[]):
+ if j and data[i, j] > 0:
+ label = N[j-1]
+ if label > k:
+ k = label
+ index = j
+ return k, index
+
+ k, kindex = max(rowmax(0, Dcols, B), colmax(0, Drows, N))
+ if k is None:
+ break
+
+ if k in B:
+ r, Br = k, kindex
+ s, Ns = colmax (Br, Drows, N)
+ if s is None:
+ raise LPError("LP problem is inconsistent")
+
+ else: # k in N
+ s, Ns = k, kindex
+ r, Br = rowmax(Ns, Dcols, B)
+ if r is None:
+ raise LPError("LP problem is dual inconsistent")
+
+ B[Br-1] = s
+ N[Ns-1] = r
+ id_rs = div(1, data[Br, Ns])
+ negid_rs = -id_rs
+ newrows = copy.deepcopy(Drows)
+
+ srows = Dcols[Ns]
+ srows.discard(Br)
+ rcols = Drows[Br]
+ rcols.discard(Ns)
+
+ for i in srows:
+ cols = newrows[i]
+ cols.update(rcols)
+ cols.discard(Ns)
+
+ for Bi in srows:
+ for Nj in newrows[Bi]:
+ d_rj = data.get((Br, Nj))
+ if d_rj is not None:
+ d_ij = data.get((Bi,Nj), 0)
+ d_ij -= d_rj*data[Bi,Ns]*id_rs
+ if d_ij==0:
+ del data[Bi,Nj]
+ else:
+ data[Bi,Nj] = d_ij
+
+ for Bi in srows:
+ data[Bi,Ns] *= id_rs
+
+ for Nj in rcols:
+ data[Br,Nj] *= negid_rs
+
+ data[Br,Ns] = id_rs
+
+ optimal = Matrix(n-1,1)
+ for i,k in enumerate(B):
+ if k>0:
+ optimal[k-1,0] = data.get((i+1,0), 0)
+ optimal_value = data.get ((0,0),0)
+ return optimal, optimal_value
+
=======================================
--- /trunk/sympycore/__init__.py Thu Nov 19 04:18:40 2009
+++ /trunk/sympycore/__init__.py Sun Jan 15 14:44:56 2012
@@ -107,6 +107,7 @@
import os
import sys
import nose
+ #mpmath = sys.modules['sympycore.arithmetic.mpmath']
d = os.path.dirname(os.path.abspath(__file__))
cmd = '%s %s %s %s' % (sys.executable, nose.core.__file__, d,
nose_args)
print >>sys.stderr, 'Running %r' % cmd
@@ -159,7 +160,9 @@
mpmath = sys.modules['sympycore.arithmetic.mpmath']
fn = os.path.join(os.path.dirname(mpmath.__file__), 'REVISION')
- l = ['mode=%s' % (mpmath.settings.MODE)]
+ l = []
+ backend = mpmath.libmp.backend.BACKEND
+ l.append ('backend=%s' % (backend))
if os.path.isfile(fn):
f = open(fn, 'r')
l.append('revision=%s' % (f.read().strip()))
@@ -169,8 +172,9 @@
print >>sys.stderr, 'mpmath version: %s (%s)' %
(mpmath.__version__, ', '.join(l))
print >>sys.stderr, 'mpmath is installed in %s' %
(os.path.dirname(mpmath.__file__))
- if mpmath.settings.MODE=='gmpy':
- gmpy = mpmath.settings.gmpy
+
+ if backend=='gmpy':
+ gmpy = mpmath.libmp.backend.gmpy
print >>sys.stderr, 'gmpy version: %s' % (gmpy.version())
print >>sys.stderr, 'gmpy is installed in %s' %
(os.path.dirname(gmpy.__file__))
=======================================
--- /trunk/sympycore/matrices/algebra.py Fri Apr 8 01:57:00 2011
+++ /trunk/sympycore/matrices/algebra.py Sun Jan 15 14:44:56 2012
@@ -603,7 +603,7 @@
return self[key, :]
elif is_integer(key):
return self[int(key)]
-
+
raise IndexError('index must be int, slice, or tuple, got %s' %
(tkey))
def _set_diagonal(self, key, value):
@@ -1030,6 +1030,7 @@
MATRIX_DICT_get_gauss_jordan_elimination_operations,
MATRIX_DICT_apply_row_operations)
from .linalg_determinant import MATRIX_DICT_determinant
+from .linalg_lp import MATRIX_DICT_LP_solve
MatrixDict.__iadd__ = MATRIX_DICT_iadd
MatrixDict.__imul__ = MATRIX_DICT_imul
@@ -1042,3 +1043,4 @@
MatrixDict.trace = MATRIX_DICT_trace
MatrixDict.get_gauss_jordan_elimination_operations =
MATRIX_DICT_get_gauss_jordan_elimination_operations
MatrixDict.apply_row_operations = MATRIX_DICT_apply_row_operations
+MatrixDict.LP_solve = MATRIX_DICT_LP_solve
=======================================
--- /trunk/sympycore/matrices/linalg.py Mon Apr 4 05:22:31 2011
+++ /trunk/sympycore/matrices/linalg.py Sun Jan 15 14:44:56 2012
@@ -693,6 +693,8 @@
return 0, row_pivot_table, pivot_table
def get_rc_map(data):
+ """ Return mapping between row indices and column indices with
existing values.
+ """
rows = {}
for i, j in data:
s = rows.get (i)
@@ -711,6 +713,7 @@
return rows
def get_rc_maps(data):
+
rows = {}
cols = {}
for i, j in data:
=======================================
--- /trunk/sympycore/matrices/polyhedra.py Sun Jan 8 14:39:40 2012
+++ /trunk/sympycore/matrices/polyhedra.py Sun Jan 15 14:44:56 2012
@@ -1,21 +1,118 @@
-
-from sympycore import Logic, SymbolicEquality, Calculus, heads
+
+#from __future__ import division
+
+import copy
+from sympycore import Logic, SymbolicEquality, Calculus, heads, Symbol
+from sympycore.arithmetic.numbers import div
from . import Matrix
class Polyhedron(object):
+ """ A basis class that provides a H-representation for general
+ polyhedra.
+
+ Consider the following polyhedron::
+
+ P(A_I, b_I, A_L, b_L) = { x in R^n | A_I*x <= b_I, A_L*x == b_L }
+
+ where ``A_I`` and ``A_L`` are matrices with ``|I|`` and ``|L|`` rows,
+ respectively, and ``I``, ``L`` are partitions of ``{0,..,m-1}`` such
that
+ ``set(I + L)==set(range(m))``, and ``*`` denotes matrix dot product.
+
+ The class ``Polyhedron`` stores inequalities and equalities in a
+ homogenized form, that is, in a m times (n+1) matrix A where the
+ first column corresponds to coefficients in ``-b_I`` and ``-b_L``,
+ and the rest of the columns correspond to coefficents in ``A_I``
+ and ``A_L``. With respect to coefficent storage, this class
+ represents the following homogeneous polyhedron::
+
+ P(A) = { x' in R^{n+1} | A[I] * x' <= 0, A[L] * x' == 0}
+
+ where ``x'`` may contain special variables, see ``labels``
+ attributes for details.
+
+ Attributes
+ ----------
+ A : Matrix
+ A constraints matrix.
+
+ I, L : list
+ Lists of row indices of ``A``. If ``i`` in ``I`` then the
+ ``i``-th row corresponds to inequality relation and if ``i`` in
+ ``L`` then to equality relation. See above for definitions.
+
+ labels : list
+ List of column labels of matrix A. Some labels have special
+ meaning, as defined below:
+
+ ``"1"`` --- the name corresponds to homogenization variable
+ and the corresponding column in A corresponds to coefficents
+ in ``-b_I`` and ``-b_L``. To relate P(A) to P(A_I, b_I, A_L,
+ b_L), ``"1"`` value should be taken to 1.
+
+ ``"SLACK#"`` --- the name corresponds to slack variable. The
+ slack variable is used in a LP problem where it is assumed
+ to be non-negative. Here ``"#"`` denotes integers used to
+ index slack variables corresponding to inequality rows.
+
+ ``"MAX"`` or ``"MIN"`` --- the labels correspond to LP
+ objective function values when used in an equality relation.
+ Minimization LP problem is automatically converted to
+ maximization LP problem.
+
+ objective_row, objective_column : {None, int}
+ Row and column indices of the objective coefficients and ``"MAX"``
+ column name, respectively.
+
+ rhs_column : {None, int}
+ Column index corresponding to RHS coefficients.
+ """
def __init__(self, *exprs):
- self.A = Matrix(0,1)
+ self.A = Matrix(0,0)
self.I = []
self.L = []
- self.names = ['*']
+ self.labels = []
+ self.objective_row = None
+ self.objective_column = None
+ self.rhs_column = None
for expr in exprs:
self.add(expr)
-
+
+ def copy(self):
+ cpy = Polyhedron()
+ cpy.A = self.A.copy()
+ for n in ['I', 'L', 'labels']:
+ getattr(cpy, n).extend(getattr(self, n))
+ for n in ['objective_row', 'objective_column']:
+ setattr(cpy, n, getattr(self, n))
+ return cpy
+
+ @property
+ def A_L(self):
+ """ Return a matrix of equalities.
+ """
+ return self.A[tuple(self.L),:]
+
+ @property
+ def A_I(self):
+ """ Return a matrix of inequalities.
+ """
+ return self.A[tuple(self.I),:]
+
def add(self, expr):
- """ Add extra contraint to polyhedron.
+ """ Add extra constraint to polyhedron.
+
+ Patameters
+ ----------
+ expr : {str, Logic}
+
+ Symbolic expression of a constraint. The expression must be
+ an equality or inequality relation of arbitrary linear
+ combination of variables.
"""
if isinstance(expr, str):
+ if not expr or expr.strip ().startswith ('#'):
+ return
if '=' in expr or '<' in expr or '>' in expr:
expr = Logic(expr)
else:
@@ -32,54 +129,127 @@
expr = lhs - rhs
self.L.append(i)
elif op in [heads.LT, heads.LE]:
- op = heads.GE
- expr = rhs - lhs
+ op = heads.LE
+ expr = lhs - rhs
self.I.append(i)
elif op in [heads.GT, heads.GE]:
- op = heads.GE
- expr = lhs - rhs
+ op = heads.LE
+ expr = rhs - lhs
self.I.append(i)
else:
raise NotImplementedError(`op`) # expected ==, <, <=, >, >=
expr = expr.to(heads.TERM_COEFF_DICT)
- rowdict = dict()
-
+ negate_row = False
if expr.head is heads.NUMBER:
- self.A[i,0] = expr.data
+ try:
+ j = self.labels.index('1')
+ except ValueError:
+ j = len(self.labels)
+ self.labels.append(Symbol('1'))
+ if self.rhs_column is None:
+ self.rhs_column = j
+ else:
+ assert self.rhs_column == j,`self.rhs_column,j`
+ self.A[i,j] = expr.data
elif expr.head is heads.TERM_COEFF:
term, coeff = expr.data
+ if term=='MIN':
+ negate_row = True
+ coeff = -coeff
+ term = Symbol('MAX')
try:
- j = self.names.index(term)
+ j = self.labels.index(term)
except ValueError:
- j = len(self.names)
- self.names.append(term)
+ j = len(self.labels)
+ self.labels.append(term)
+ if term=='MAX':
+ assert self.objective_row is None, `self.objective_row,
term`
+ self.objective_row = i
+ self.objective_column = j
self.A[i,j] = coeff
elif expr.head is heads.TERM_COEFF_DICT:
for term, coeff in expr.data.iteritems():
if term.head == heads.NUMBER:
assert term.data==1,`term.pair` # expected
Calculus('1')
- self.A[i,0] = coeff
+ try:
+ j = self.labels.index('1')
+ except ValueError:
+ j = len(self.labels)
+ self.labels.append(Symbol('1'))
+ if self.rhs_column is None:
+ self.rhs_column = j
+ else:
+ assert self.rhs_column == j,`self.rhs_column,j`
+ self.A[i,j] = coeff
else:
+ if term=='MIN':
+ negate_row = True
+ coeff = -coeff
+ term = Symbol('MAX')
try:
- j = self.names.index(term)
+ j = self.labels.index(term)
except ValueError:
- j = len(self.names)
- self.names.append(term)
+ j = len(self.labels)
+ self.labels.append(term)
+ if term=='MAX':
+ assert self.objective_row is None,
`self.objective_row, term`
+ self.objective_row = i
+ self.objective_column = j
self.A[i,j] = coeff
else:
try:
- j = self.names.index(expr)
+ j = self.labels.index(expr)
except ValueError:
- j = len(self.names)
- self.names.append(expr)
- self.A[i,j] = 1
-
- self.A = self.A.resize(i+1, len(self.names))
+ j = len(self.labels)
+ self.labels.append(expr)
+ self.A[i,j] = 1
+
+ self.A = self.A.resize(i+1, len (self.labels))
+
+ def get_LP(self):
+ """Construct LP dictionary matrix from polyhedra constraints.
+
+ Parameters
+ ----------
+ None
+
+ Returns
+ -------
+ variables : list
+ List of variables.
+ D : Matrix
+ LP dictionary matrix.
+
+ See also
+ --------
+ MatrixDict.LP_solve
+ """
+ assert len (self.L)==1,`self.L, self.objective_row`
+ colindices = []
+
+ variables = []
+ for name, i in sorted (zip(map (str, self.labels), range (len
(self.labels)))):
+ if name=='1':
+ pass
+ elif name=='MAX':
+ pass
+ else:
+ colindices.append(i)
+ variables.append(name)
+ colindices = tuple(colindices)
+ rowindices = tuple(self.I)
+ D = Matrix(self.A.rows, len (colindices)+1)
+ D[1:,1:] = -self.A[rowindices,colindices]
+ D[0,1:] = self.A[self.objective_row, colindices] /
(-self.A[self.objective_row, self.objective_column])
+ if self.rhs_column is not None:
+ D[:,0] = -self.A[:,self.rhs_column]
+
+ return variables, D
def show(self):
- print 'A:',', '.join(map(str, self.names))
+ print 'A:',', '.join(map(str, self.labels))
print self.A
print 'L:',self.L
print 'I:',self.I
=======================================
--- /trunk/sympycore/matrices/tests/test_polyhedra.py Sun Jan 8 14:39:40
2012
+++ /trunk/sympycore/matrices/tests/test_polyhedra.py Sun Jan 15 14:44:56
2012
@@ -1,8 +1,49 @@
-
-from sympycore import Polyhedron
+
+import time
+from sympycore import Polyhedron, Matrix, mpq
+
+
+def test_LP_solve_Chateau_ETH_Production():
+ constraints = '''\
+3*x1+4*x2+2*x3==MAX
+2*x1<=4
+x1+2*x3<=8
+3*x2+x3<=6
+'''
+ p = Polyhedron(*constraints.split('\n'))
+ #p.show()
+ names, D = p.get_LP()
+ Dwork = D[:]
+ xopt, vopt = Dwork.LP_solve(overwrite=True)
+ assert (D[0,1:] * xopt)[0,0]==vopt==16,`D[0,1:] * xopt,vopt`
+
+def test_LP_solve_cyclic():
+ constraints = '''
+x1-2*x2+x3==MAX
+2*x1-x2+x3<=0
+3*x1+x2+x3<=0
+-5*x1+3*x2-2*x3<=0
+'''
+ p = Polyhedron(*constraints.split('\n'))
+ #p.show()
+ names, D = p.get_LP()
+ Dwork = D[:]
+ xopt, vopt = Dwork.LP_solve(overwrite=True)
+ assert (D[0,1:] * xopt)[0,0]==vopt==0,`D[0,1:] * xopt,vopt`
+
+def test_LP_solve_Karmarkar_example ():
+ p = Polyhedron('x1+x2==MAX')
+ for i in range (11):
+ p1 = mpq((i,10))*1
+ p.add('2*%s*x1+x2<=%s+1' % (p1,p1**2))
+ names, D = p.get_LP()
+ Dwork = D[:]
+ xopt, vopt = Dwork.LP_solve(overwrite=True)
+ assert (D[0,1:] * xopt)[0,0]==vopt==mpq ((5,4)),`D[0,1:] * xopt,vopt`
def test_Av98a():
constraints = '''\
+#x1+x2+x3==1
1+x1>=0
1+x2>=0
1-x1>=0
@@ -12,6 +53,30 @@
1+x1+x3>=0
1+x2+x3>=0'''
p = Polyhedron(*constraints.split('\n'))
- p.show ()
+ #p.show()
+
+ return
+
vertices = [(1,1,0),(-1,1,0),(1,-1,0),(-1,-1,0),(0,0,-1)]
rays = [(0,0,1)]
+
+ print 'A_L'
+ print p.A_L
+
+ xd, dep, indep = p.A_L.solve_null(labels=p.names)
+ Aker = Matrix([xd[l].tolist()[0] for l in p.names])
+ print 'Aker:'
+ print Aker
+
+ print 'A * Aker:'
+ print p.A * Aker
+
+ xd, dep, indep = p.A_I[(0,1),:].solve_null(labels=p.names)
+ Aker = Matrix([xd[l].tolist()[0] for l in p.names])
+ print 'Aker:'
+ print Aker
+ print 'A * Aker:'
+ print p.A * Aker
+
+ print p.A[(0,1,2,4),:].I*2
+ print p.A[(0,1,2,4),:]*p.A[(0,1,2,4),:].I