denis-bz/0-lpgen.md

## 0-lpgen.md

      
    Raw
  

              0-lpgen.md
            
          
    Test generator for sparse linear programming: m+n constraints, mn variables, 2mn non-zeros

Keywords: linear-programming, sparse, scipy, test-case-generator
Purpose

lpgen_2d generates test cases for linear programming, with parameters m and n:
m+n constraints, mn variables, 2mn non-zeros. An example:
from scipy.optimize import linprog

A, b, c = lpgen_2d( m, n, sparse=True, dtype=np.float32, seed=0, verbose=1 )
    # e.g.
    # lpgen_2d: m 1000  n 2000  A (3000, 2000000) 4e+07 bytes sparse float32  seed 0

res = linprog( c, A_ub=A, b_ub=b, method="interior-point" ... )
    # takes ~ 90 sec real time, 120 cpu time (sum of cores) on my 2.7 GHz iMac

    Sec  Objective    Primal, Dual Feasibility   Gap   Step   Path Parameter
      0: -1.00036e+07 1          1          1          nan        1         
     21: -8.53131e+06 0.0013     0.0013     0.0013     0.999      0.0013    
     41: -75664.5     0.000467   0.000467   0.000467   0.707      0.000467  
     ...
     122: -20000      1.12e-11   1.12e-11   1.12e-11   1          1.12e-11  

Notes

Lpgen_2d with integer n/m generates
assignment problems,
which have many vertices, so make pretty good test problems for LP.
The simplex method finds integer solutions;
interior-point generally does not, but is faster by a huge factor.
For example, for m=50 n=100 ip takes 0.2 seconds (real), simplex 13 seconds.
Only linprog( method="interior-point" ) handles sparse matrices;
the default method="simplex" hits bugs (in scipy 1.2.1, April 2019).
All LP solvers have many many options ...
Sparse matrices are only ~ 1/4 smaller with float32 than float64 --
indices are 64-bit ints.
Multicore ? I see cpu and wall clock times e.g. 145 sec, 90 sec
(time.clock sum all cores, time.time).
The source has a "TODO: if umfpack is installed ..."
See also

Scipy linprog interior-point

_linprog_ip.py  -- 1100 lines, well commented

Assignment problem in
Knuth, Stanford Graphbase
pages 26-30, 77-87.

Comments welcome, more test cases welcome.

cheers

-- denis 27 April 2019

  
## lpgen.py
#!/usr/bin/env python
""" lpgen_2d LP test: A m+n constraints, m*n vars, density 2/(m+n) """

from __future__ import division, print_function
import numpy as np
from scipy import sparse as sp

__version__ = "2019-04-23 april  denis-bz-py"  # sparse

#...............................................................................
def lpgen_2d( m, n, sparse=True, dtype=np.float32, seed=0, cint=10, verbose=1 ):
    """ -> A b c LP test: A (m+n constraints, m*n vars), density 2/(m+n)
        row sums == n/m, col sums == 1
        c random integers if cint / random exponential
        sparse: linprog method="interior-point" only
    """
        # test-linprog-2d.py m=1000 sparse: float32 float64 40m 56m, both 90 sec
    m = int(m)
    n = int(n)
    random = seed if isinstance(seed, np.random.RandomState) \
            else np.random.RandomState( seed=seed )
    if cint > 0:
        c = - random.randint( 0, cint+1, size=(m,n) ).astype(dtype)
    else:
        c = - random.exponential( size=(m,n) ).astype(dtype)
    zeros = np.zeros if not sparse \
        else sp.lil_matrix
    Arow = zeros(( m, m*n ), dtype=dtype )
    Acol = zeros(( n, m*n ), dtype=dtype )
    brow = np.zeros( m, dtype=dtype )
    bcol = np.zeros( n, dtype=dtype )
    for j in range(m):
        j1 = j + 1
        Arow[ j, j*n : j1*n ] = 1
        brow[j] = n / m
    for j in range(n):
        Acol[ j, j::n ] = 1
        bcol[j] = 1
    if not sparse:
        A = np.vstack(( Arow, Acol ))
        nbytes = A.nbytes
    else:
        A = sp.vstack(( Arow, Acol )).tocsc()
        nbytes = A.data.nbytes + A.indices.nbytes + A.indptr.nbytes
            # float32 ~ 8+4 bytes per val, float64 ~ 8+8
    b = np.hstack(( brow, bcol ))
    c = c.reshape(-1)
    if verbose:
        sparsedense = "sparse" if sparse else "dense"
        print( "lpgen_2d: m %d  n %d  A %s %.3g bytes %s %s  seed %d  -c %s ..." % (
                m, n, A.shape, nbytes, sparsedense, dtype.__name__, seed, - c[:10] ))
    return A, b, c


def _ints( x ):
    return np.asanyarray(x) .round().astype(int)

#...............................................................................
if __name__ == "__main__":
    import sys

    np.set_printoptions( threshold=20, edgeitems=10, linewidth=140,
            formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g

    m = 2
    n = 4
    sparse = True
    seed = 0

    # to change these params, run this.py  a=1  b=None  c='expr' ...
    # in sh or ipython --
    for arg in sys.argv[1:]:
        exec( arg )

    A, b, c = lpgen_2d( m, n, sparse=sparse, seed=seed, verbose=2 )


## test-linprog-2d-m1000.log
# from: test-linprog-2d.py m=1000 sparse=1 dtype=np.float32 methods=[0]
# 23 Apr 2019 14:35  ~bz/py/opt/lp/lpgen  Denis-iMac 10.10.5
versions: numpy 1.16.3  scipy 1.2.1  python 2.7.15  mac 10.10.5


================================================================================
params: m 1000  c 2  bounds None  eq False  tol 0.0001  sparse  dtype float32  seed 0
test-linprog-2d.py run in ~bz/py/opt/lp/lpgen  2019-04-23 Apr 14:35
versions: numpy 1.16.3  scipy 1.2.1   python 2.7.15   mac 10.10.5
lpgen_2d: m 1000  n 2000  A (3000, 2000000) 4e+07 bytes sparse float32  seed 0  -c [5 0 3 3 7 9 3 5 2 4] ...
Sec  Objective    Primal, Dual Feasibility   Gap   Step   Path Parameter
  0: -1.00036e+07 1          1          1          nan        1
 21: -8.53131e+06 0.0013     0.0013     0.0013     0.999      0.0013
 41: -75664.5     0.000467   0.000467   0.000467   0.707      0.000467
 61: -28240.4     0.000124   0.000124   0.000124   0.743      0.000124
 81: -21904.6     3.09e-05   3.09e-05   3.09e-05   0.827      3.09e-05
102: -20013.3     2.23e-07   2.23e-07   2.23e-07   0.994      2.23e-07
122: -20000       1.12e-11   1.12e-11   1.12e-11   1          1.12e-11
Optimization terminated successfully.
         Current function value: -20000.000670
         Iterations: 6

test_linprog interior-point: f -20000  niter 6  A (3000, 2000000)  bounds None  eq False
test_linprog: 135298 x are in .1 .. .99  min av max 0.01 0.0115 0.439
test_linprog: b - Ax min av max -1.12e-07 -7.46e-08 -5.57e-08
89.0 145.5 sec  interior-point  m 1000  n 2000  sparse: niter 6
b - A . x.round  quantiles: 1 1 1 2 2  counts: [2000 1000]
       92.05 real       129.76 user        19.34 sys

## test-linprog-2d.py
#!/usr/bin/env python2
""" test-linprog-2d.py: lpgen_2d, test_linprog """

from __future__ import division, print_function
import numpy as np
from scipy.optimize import linprog  # $scopt/_linprog.py linprog.doc
    # https://docs.scipy.org/doc/scipy/reference/optimize.linprog-interior-point.html

from lpgen import lpgen_2d  # -> A b c: m+n constraints, m*n vars

__version__ = "2019-04-25 april  denis-bz-py"  # sparse

#...............................................................................
def test_linprog( A, b, c, method="interior-point", bounds=None, eq=False, **kw ):
    """ linprog( c, A, b ) eq / ub  ->  linprog """
    if eq:
        res = linprog( c, A_eq=A, b_eq=b, method=method, bounds=bounds, **kw )
    else:
        res = linprog( c, A_ub=A, b_ub=b, method=method, bounds=bounds, **kw )
    if res.status != 0:
        print( res.message )
        return res
    x = res.x
    print( "\ntest_linprog %s: f %g  niter %d  A %s  bounds %s  eq %s " % (
            method, res.fun, res.nit, A.shape, bounds, eq ))

        # simplex all 0-1, interior some in between, b - Ax ?
    nz = x[ (.01 < x) & (x < .99) ]
    if len(nz) > 0:
        print( "test_linprog: %d x are in .1 .. .99  %s" % (
                len(nz), _minavmax( nz )))
    b_Ax = b - A.dot( x )  # >= 0 ?
    print( "test_linprog: b - Ax", _minavmax( b_Ax ))

    return res  # an OptimizeResult


def _minavmax( x ):
    return "min av max %.3g %.3g %.3g" % (x.min(), x.mean(), x.max())


#...............................................................................
if __name__ == "__main__":  # test_linprog_2d( m, n ) for m n seed ...
    import sys
    import z_numpyutil as nu

    np.set_printoptions( threshold=20, edgeitems=20, linewidth=140,
            formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g
    print( "\n" + 80 * "=" )

        # lpgen: A (m+n constraints, m*n vars), density 2/(m+n)
    m = 50  # row sums == c, col sums == 1
    c = 2  # n = c * m
    sparse = True  # bug _linprog_simplex 566 A[is_negative_constraint] *= -1
    dtype = np.float32  # m 1000 32 64: 40 56 mbytes, both ~ 120 sec
    seed = 0

        # linprog --
    methodsdict = { 0: "interior-point", 1: "simplex" }
    methods = [0]
    bounds = None  # == np.r_[ 0, np.inf ] ?
    eq = False
        # options --
    disp = True  # patch $scipy/optimize/_linprog_ip.py clock, alpha NaN not "-"
    maxiter = 1e5
    tol = 1e-4

    # to change these params, run this.py  a=1  b=None  c='expr' ...
    # in sh or ipython --
    for arg in sys.argv[1:]:
        exec( arg )

    np.random.seed(seed)
    sparsedense = "sparse" if sparse else "dense"
    print( "params: m %d  c %g  bounds %s  eq %s  tol %g  %s  dtype %s  seed %s " % (
            m, c, bounds, eq, tol, sparsedense, dtype.__name__, seed ))
    options = dict( disp=disp, maxiter=int(maxiter), sparse=sparse, tol=tol )

    print( nu.From( __file__ ))  # file pwd date
    print( nu.versions() )

    n = int( m * c )
    A, b, c = lpgen_2d( m, n, sparse=sparse, dtype=dtype, seed=seed )

    for jmethod in methods:
        method = methodsdict[jmethod]
        if method == "simplex" and sparse:
            print( "? linprog( sparse A, simplex ) may be buggy" )
            continue
        nu.ptime()
#...............................................................................
        res = test_linprog( A, b, c, method=method, bounds=bounds, eq=eq,
                options=options )

        nu.ptime( "%s  m %d  n %d  %s: niter %d " % (
                method, m, n, sparsedense, res.nit ))
        if method == "interior-point":  # b - A . x.round ?
            b_Ax = b - A.dot( res.x.round() )
            q = nu.quantiles( b_Ax, q=5 )
            ints = b_Ax.round().astype(int)
            counts = np.bincount( ints - ints.min() )  # m 50 c 2: [ 1 28 95 26]
            print( "b - A . x.round  %s  counts: %s " % (q, counts) )
	#!/usr/bin/env python
	""" lpgen_2d LP test: A m+n constraints, m*n vars, density 2/(m+n) """

	from __future__ import division, print_function
	import numpy as np
	from scipy import sparse as sp

	__version__ = "2019-04-23 april denis-bz-py" # sparse

	#...............................................................................
	def lpgen_2d( m, n, sparse=True, dtype=np.float32, seed=0, cint=10, verbose=1 ):
	""" -> A b c LP test: A (m+n constraints, m*n vars), density 2/(m+n)
	row sums == n/m, col sums == 1
	c random integers if cint / random exponential
	sparse: linprog method="interior-point" only
	"""
	# test-linprog-2d.py m=1000 sparse: float32 float64 40m 56m, both 90 sec
	m = int(m)
	n = int(n)
	random = seed if isinstance(seed, np.random.RandomState) \
	else np.random.RandomState( seed=seed )
	if cint > 0:
	c = - random.randint( 0, cint+1, size=(m,n) ).astype(dtype)
	else:
	c = - random.exponential( size=(m,n) ).astype(dtype)
	zeros = np.zeros if not sparse \
	else sp.lil_matrix
	Arow = zeros(( m, m*n ), dtype=dtype )
	Acol = zeros(( n, m*n ), dtype=dtype )
	brow = np.zeros( m, dtype=dtype )
	bcol = np.zeros( n, dtype=dtype )
	for j in range(m):
	j1 = j + 1
	Arow[ j, jn : j1n ] = 1
	brow[j] = n / m
	for j in range(n):
	Acol[ j, j::n ] = 1
	bcol[j] = 1
	if not sparse:
	A = np.vstack(( Arow, Acol ))
	nbytes = A.nbytes
	else:
	A = sp.vstack(( Arow, Acol )).tocsc()
	nbytes = A.data.nbytes + A.indices.nbytes + A.indptr.nbytes
	# float32 ~ 8+4 bytes per val, float64 ~ 8+8
	b = np.hstack(( brow, bcol ))
	c = c.reshape(-1)
	if verbose:
	sparsedense = "sparse" if sparse else "dense"
	print( "lpgen_2d: m %d n %d A %s %.3g bytes %s %s seed %d -c %s ..." % (
	m, n, A.shape, nbytes, sparsedense, dtype.__name__, seed, - c[:10] ))
	return A, b, c


	def _ints( x ):
	return np.asanyarray(x) .round().astype(int)

	#...............................................................................
	if __name__ == "__main__":
	import sys

	np.set_printoptions( threshold=20, edgeitems=10, linewidth=140,
	formatter = dict( float = lambda x: "%.2g" % x )) # float arrays %.2g

	m = 2
	n = 4
	sparse = True
	seed = 0

	# to change these params, run this.py a=1 b=None c='expr' ...
	# in sh or ipython --
	for arg in sys.argv[1:]:
	exec( arg )

	A, b, c = lpgen_2d( m, n, sparse=sparse, seed=seed, verbose=2 )
	# from: test-linprog-2d.py m=1000 sparse=1 dtype=np.float32 methods=[0]
	# 23 Apr 2019 14:35 ~bz/py/opt/lp/lpgen Denis-iMac 10.10.5
	versions: numpy 1.16.3 scipy 1.2.1 python 2.7.15 mac 10.10.5


	================================================================================
	params: m 1000 c 2 bounds None eq False tol 0.0001 sparse dtype float32 seed 0
	test-linprog-2d.py run in ~bz/py/opt/lp/lpgen 2019-04-23 Apr 14:35
	versions: numpy 1.16.3 scipy 1.2.1 python 2.7.15 mac 10.10.5
	lpgen_2d: m 1000 n 2000 A (3000, 2000000) 4e+07 bytes sparse float32 seed 0 -c [5 0 3 3 7 9 3 5 2 4] ...
	Sec Objective Primal, Dual Feasibility Gap Step Path Parameter
	0: -1.00036e+07 1 1 1 nan 1
	21: -8.53131e+06 0.0013 0.0013 0.0013 0.999 0.0013
	41: -75664.5 0.000467 0.000467 0.000467 0.707 0.000467
	61: -28240.4 0.000124 0.000124 0.000124 0.743 0.000124
	81: -21904.6 3.09e-05 3.09e-05 3.09e-05 0.827 3.09e-05
	102: -20013.3 2.23e-07 2.23e-07 2.23e-07 0.994 2.23e-07
	122: -20000 1.12e-11 1.12e-11 1.12e-11 1 1.12e-11
	Optimization terminated successfully.
	Current function value: -20000.000670
	Iterations: 6

	test_linprog interior-point: f -20000 niter 6 A (3000, 2000000) bounds None eq False
	test_linprog: 135298 x are in .1 .. .99 min av max 0.01 0.0115 0.439
	test_linprog: b - Ax min av max -1.12e-07 -7.46e-08 -5.57e-08
	89.0 145.5 sec interior-point m 1000 n 2000 sparse: niter 6
	b - A . x.round quantiles: 1 1 1 2 2 counts: [2000 1000]
	92.05 real 129.76 user 19.34 sys
	#!/usr/bin/env python2
	""" test-linprog-2d.py: lpgen_2d, test_linprog """

	from __future__ import division, print_function
	import numpy as np
	from scipy.optimize import linprog # $scopt/_linprog.py linprog.doc
	# https://docs.scipy.org/doc/scipy/reference/optimize.linprog-interior-point.html

	from lpgen import lpgen_2d # -> A b c: m+n constraints, m*n vars

	__version__ = "2019-04-25 april denis-bz-py" # sparse

	#...............................................................................
	def test_linprog( A, b, c, method="interior-point", bounds=None, eq=False, **kw ):
	""" linprog( c, A, b ) eq / ub -> linprog """
	if eq:
	res = linprog( c, A_eq=A, b_eq=b, method=method, bounds=bounds, **kw )
	else:
	res = linprog( c, A_ub=A, b_ub=b, method=method, bounds=bounds, **kw )
	if res.status != 0:
	print( res.message )
	return res
	x = res.x
	print( "\ntest_linprog %s: f %g niter %d A %s bounds %s eq %s " % (
	method, res.fun, res.nit, A.shape, bounds, eq ))

	# simplex all 0-1, interior some in between, b - Ax ?
	nz = x[ (.01 < x) & (x < .99) ]
	if len(nz) > 0:
	print( "test_linprog: %d x are in .1 .. .99 %s" % (
	len(nz), _minavmax( nz )))
	b_Ax = b - A.dot( x ) # >= 0 ?
	print( "test_linprog: b - Ax", _minavmax( b_Ax ))

	return res # an OptimizeResult


	def _minavmax( x ):
	return "min av max %.3g %.3g %.3g" % (x.min(), x.mean(), x.max())


	#...............................................................................
	if __name__ == "__main__": # test_linprog_2d( m, n ) for m n seed ...
	import sys
	import z_numpyutil as nu

	np.set_printoptions( threshold=20, edgeitems=20, linewidth=140,
	formatter = dict( float = lambda x: "%.2g" % x )) # float arrays %.2g
	print( "\n" + 80 * "=" )

	# lpgen: A (m+n constraints, m*n vars), density 2/(m+n)
	m = 50 # row sums == c, col sums == 1
	c = 2 # n = c * m
	sparse = True # bug _linprog_simplex 566 A[is_negative_constraint] *= -1
	dtype = np.float32 # m 1000 32 64: 40 56 mbytes, both ~ 120 sec
	seed = 0

	# linprog --
	methodsdict = { 0: "interior-point", 1: "simplex" }
	methods = [0]
	bounds = None # == np.r_[ 0, np.inf ] ?
	eq = False
	# options --
	disp = True # patch $scipy/optimize/_linprog_ip.py clock, alpha NaN not "-"
	maxiter = 1e5
	tol = 1e-4

	# to change these params, run this.py a=1 b=None c='expr' ...
	# in sh or ipython --
	for arg in sys.argv[1:]:
	exec( arg )

	np.random.seed(seed)
	sparsedense = "sparse" if sparse else "dense"
	print( "params: m %d c %g bounds %s eq %s tol %g %s dtype %s seed %s " % (
	m, c, bounds, eq, tol, sparsedense, dtype.__name__, seed ))
	options = dict( disp=disp, maxiter=int(maxiter), sparse=sparse, tol=tol )

	print( nu.From( __file__ )) # file pwd date
	print( nu.versions() )

	n = int( m * c )
	A, b, c = lpgen_2d( m, n, sparse=sparse, dtype=dtype, seed=seed )

	for jmethod in methods:
	method = methodsdict[jmethod]
	if method == "simplex" and sparse:
	print( "? linprog( sparse A, simplex ) may be buggy" )
	continue
	nu.ptime()
	#...............................................................................
	res = test_linprog( A, b, c, method=method, bounds=bounds, eq=eq,
	options=options )

	nu.ptime( "%s m %d n %d %s: niter %d " % (
	method, m, n, sparsedense, res.nit ))
	if method == "interior-point": # b - A . x.round ?
	b_Ax = b - A.dot( res.x.round() )
	q = nu.quantiles( b_Ax, q=5 )
	ints = b_Ax.round().astype(int)
	counts = np.bincount( ints - ints.min() ) # m 50 c 2: [ 1 28 95 26]
	print( "b - A . x.round %s counts: %s " % (q, counts) )