Skip to content

Instantly share code, notes, and snippets.

@denis-bz

denis-bz/0-lpgen.md

Last active May 6, 2021 19:31
Show Gist options
  • Save denis-bz/8647461 to your computer and use it in GitHub Desktop.
Save denis-bz/8647461 to your computer and use it in GitHub Desktop.
Download ZIP
Test generator for sparse linear programming: m+n constraints, mn variables, 2mn non-zeros 26 Apr 2019
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

Raw
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 )
Raw
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
Raw
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) )
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment