login  home  contents  what's new  discussion  bug reports help  links  subscribe  changes  refresh  edit

)abbrev package LINPROG LinearProgramming
++ Author: Kurt Pagani
++ Date Created: Sun Jun 03 02:05:55 CEST 2018
++ References:
++ Description:
++ TODO: description, more in test_linprog, eps, tol
++ IDEA: linSolve standard method, try LU from jet
++ LINPROG
++
LinearProgramming() : Exports == Implementation where
R ==> Float
PI ==> PositiveInteger
NNI ==> NonNegativeInteger
MATR ==> Matrix R
VECR ==> Vector R
LI ==> List Integer
RESR ==> Record(x:VECR,f:R,itn:NNI)
RESRS ==> Record(X:VECR,objVal:R,Basis:LI)
RLUP ==> Record(L:MATR,U:MATR,P:MATR)
RXP ==> Record(X:VECR,P:MATR)
Exports ==  with
revisedSimplex : (MATR,VECR,VECR,LI) -> RESRS
revisedSimplexLU : (MATR,VECR,VECR,LI) -> RESRS
luDecompose : MATR -> RLUP
luSolve : (MATR,VECR) -> RXP
Implementation ==  MATR add
eps:R:=0.00001
zeroTol:R:= 0.00001
revisedSimplex(A:MATR,b:VECR,c:VECR,B:LI):RESRS ==
numRows:PI := #b ::PI
numVars:PI := #c ::PI
--zeroTol:R:= 0.00001
N:LI := [i for i in 1..numVars | not member?(i,B)]
--N:LI:=setDifference([i for i in 1..numVars],B)
RA:LI:=[i for i in 1..nrows(A)]
AB:MATR:=A(RA,B)
AN:MATR:=A(RA,N)
numNonbasicVars:=#N
cN := vector [c.i for i in N]
cB := vector [c.i for i in B]
IAB:Union(MATR,"failed"):=inverse(AB)
if IAB case MATR then
bBAR:VECR := IAB*b
else
error "AB not invertible"
objVal := dot(cB,bBAR)
negRedCost := -1.0
while negRedCost < 0.0 repeat
IAB:Union(MATR,"failed"):=inverse(AB)
if IAB case MATR then
wN:VECR := cN - (cB*IAB)*AN
else
error "AB not invertible"
--
negRedCostIdx:NNI := 0
negRedCost := 0.0
--
for k in 1..numNonbasicVars repeat
if wN.k < -zeroTol then
negRedCost := wN.k
--record the index
negRedCostIdx := k
break
--
if negRedCost < 0.0 then
IAB:Union(MATR,"failed"):=inverse(AB)
if IAB case MATR then
aBAR:VECR := IAB*column(AN,negRedCostIdx)
else
error "AB not invertible"
minRatioVal := 1.0/zeroTol --inf
minRatioIdx := 0
--
for k in 1..numRows repeat
if aBAR.k > zeroTol then
if bBAR.k / aBAR.k < minRatioVal then
minRatioVal := bBAR.k / aBAR.k
--record the index
minRatioIdx:NNI := k
--
bBAR := bBAR - minRatioVal * aBAR
bBAR.minRatioIdx := minRatioVal
--
tmpIdx := B.minRatioIdx
B.minRatioIdx := N.negRedCostIdx
N.negRedCostIdx := tmpIdx
--
AB:MATR:=A(RA,B)
AN:MATR:=A(RA,N)
cN := vector [c.i for i in N]
cB := vector [c.i for i in B]
X:VECR:=zero(numVars)
for i in 1..numRows repeat
X.(B.i) := bBAR.i
objVal := dot(cB,bBAR)
return [X,objVal,B]
luDecompose(A:MATR):RLUP ==
not square?(A)$MATR => error "Input matrix is not square." (numRows, numVars) := (nrows A, ncols A) U:MATR:=copy A LK:MATR:=zero(numRows,numRows)$MATR
I:MATR:=diagonalMatrix(vector [1 for i in 1..numRows])
P:MATR:=copy I
LINV:=copy I
RA:LI:=[i for i in 1..numRows]
for i in 1..numVars repeat
perm:LI:=[m for m in 1..numRows]   -- don't use RA or i as index!!
pivotEl := 0$R for j in i..numRows repeat if abs(U(j,i)) > eps then tmp := perm.i perm.i := perm.j perm.j := tmp U:=I(perm,RA)*U P:=I(perm,RA)*P pivotEl := U(i, i) break abs(pivotEl) <= eps => error("Matrix is singular") LK := I(perm,RA) * LK LK(i,i) := 1$R
LINVK := copy I
--
for k in i+1..numRows repeat
if U(k, i) ~= 0$R then multiplier := U(k, i)/pivotEl U(k,RA) := U(k,RA) - multiplier * U(i,RA) LK(k,i) := multiplier LINVK(k, i) := - multiplier -- LINV := LINVK*LINV return [LK,U,P] luSolve(A:MATR,b:VECR):RXP == not square?(A)$MATR => error "Input matrix is not square."
(numRows, numVars) := (nrows A, ncols A)
r:RLUP:=luDecompose(A)
U:=copy r.U
L:=copy r.L
P:=copy r.P
Y:=vector [0$R for m in 1..numRows] X:=vector [0$R for m in 1..numRows]
rhs:VECR := P*b
for i in 1..numVars repeat
Y.i := rhs.i
rhs := rhs - Y.i * column(L,i)
for i in numVars..1 by -1 repeat
X.i := Y.i / U(i,i)
Y := Y - X.i * column(U,i)
return [X,P]
revisedSimplexLU(A:MATR,b:VECR,c:VECR,B:LI):RESRS ==
numRows:PI := #b ::PI
numVars:PI := #c ::PI
--zeroTol:R:= 0.00001
N:LI := [i for i in 1..numVars | not member?(i,B)]
--N:LI:=setDifference([i for i in 1..numVars],B)
RA:LI:=[i for i in 1..nrows(A)]
AB:MATR:=A(RA,B)
AN:MATR:=A(RA,N)
numNonbasicVars:=#N
cN := vector [c.i for i in N]
cB := vector [c.i for i in B]
--
bBAR:VECR := luSolve(AB,b).X ---- luSolve
objVal := dot(cB,bBAR)
negRedCost := -1.0
while negRedCost < 0.0 repeat
u:VECR:= luSolve(transpose(AB),cB).X
wN:VECR := cN - u*AN
--
negRedCostIdx:NNI := 0
negRedCost := 0.0
--
for k in 1..numNonbasicVars repeat
if wN.k < -zeroTol then
negRedCost := wN.k
--record the index
negRedCostIdx := k
break
--
if negRedCost < 0.0 then
aBAR:VECR := luSolve(AB,column(AN,negRedCostIdx)).X
minRatioVal := 1.0/zeroTol --inf
minRatioIdx := 0
--
for k in 1..numRows repeat
if aBAR.k > zeroTol then
if bBAR.k / aBAR.k < minRatioVal then
minRatioVal := bBAR.k / aBAR.k
--record the index
minRatioIdx:NNI := k
--
bBAR := bBAR - minRatioVal * aBAR
bBAR.minRatioIdx := minRatioVal
--
tmpIdx := B.minRatioIdx
B.minRatioIdx := N.negRedCostIdx
N.negRedCostIdx := tmpIdx
--
AB:MATR:=A(RA,B)
AN:MATR:=A(RA,N)
cN := vector [c.i for i in N]
cB := vector [c.i for i in B]
X:VECR:=zero(numVars)
for i in 1..numRows repeat
X.(B.i) := bBAR.i
objVal := dot(cB,bBAR)
return [X,objVal,B]
   Compiling FriCAS source code from file
using old system compiler.
LINPROG abbreviates package LinearProgramming
------------------------------------------------------------------------
initializing NRLIB LINPROG for LinearProgramming
compiling into NRLIB LINPROG
compiling exported revisedSimplex : (Matrix Float,Vector Float,Vector Float,List Integer) -> Record(X: Vector Float,objVal: Float,Basis: List Integer)
Local variable minRatioIdx type redefined: (NonNegativeInteger) to (Float)
Time: 0.11 SEC.
compiling exported luDecompose : Matrix Float -> Record(L: Matrix Float,U: Matrix Float,P: Matrix Float)
Time: 0.04 SEC.
compiling exported luSolve : (Matrix Float,Vector Float) -> Record(X: Vector Float,P: Matrix Float)
Time: 0.02 SEC.
compiling exported revisedSimplexLU : (Matrix Float,Vector Float,Vector Float,List Integer) -> Record(X: Vector Float,objVal: Float,Basis: List Integer)
Local variable minRatioIdx type redefined: (NonNegativeInteger) to (Float)
Time: 0.03 SEC.
(time taken in buildFunctor:  0)
;;;     ***       |LinearProgramming| REDEFINED
;;;     ***       |LinearProgramming| REDEFINED
Time: 0.01 SEC.
Warnings:
[1] revisedSimplex:  negRedCostIdx has no value
[2] revisedSimplex:  minRatioIdx has no value
[3] revisedSimplex:  bBAR has no value
[4] revisedSimplex:  X has no value
[5] luDecompose:  U has no value
[6] luDecompose:  LK has no value
[7] luDecompose:  P has no value
[8] luSolve:  U has no value
[9] luSolve:  L has no value
[10] luSolve:  P has no value
[11] revisedSimplexLU:  negRedCostIdx has no value
[12] revisedSimplexLU:  minRatioIdx has no value
[13] revisedSimplexLU:  bBAR has no value
[14] revisedSimplexLU:  X has no value
Cumulative Statistics for Constructor LinearProgramming
Time: 0.21 seconds
--------------non extending category----------------------
.. LinearProgramming of cat
(CATEGORY |package|
(SIGNATURE |revisedSimplex|
((|Record| (|:| X (|Vector| (|Float|))) (|:| |objVal| (|Float|))
(|:| |Basis| (|List| (|Integer|))))
(|Matrix| (|Float|)) (|Vector| (|Float|)) (|Vector| (|Float|))
(|List| (|Integer|))))
(SIGNATURE |revisedSimplexLU|
((|Record| (|:| X (|Vector| (|Float|))) (|:| |objVal| (|Float|))
(|:| |Basis| (|List| (|Integer|))))
(|Matrix| (|Float|)) (|Vector| (|Float|)) (|Vector| (|Float|))
(|List| (|Integer|))))
(SIGNATURE |luDecompose|
((|Record| (|:| L (|Matrix| (|Float|))) (|:| U (|Matrix| (|Float|)))
(|:| P (|Matrix| (|Float|))))
(|Matrix| (|Float|))))
(SIGNATURE |luSolve|
((|Record| (|:| X (|Vector| (|Float|))) (|:| P (|Matrix| (|Float|))))
(|Matrix| (|Float|)) (|Vector| (|Float|)))))   has no
(|MatrixCategory| (|Float|) (|Vector| (|Float|)) (|Vector| (|Float|)))    finalizing NRLIB LINPROG
Processing LinearProgramming for Browser database:
--------constructor---------
--->-->LinearProgramming((revisedSimplex ((Record (: X (Vector (Float))) (: objVal (Float)) (: Basis (List (Integer)))) (Matrix (Float)) (Vector (Float)) (Vector (Float)) (List (Integer))))): Not documented!!!!
--->-->LinearProgramming((revisedSimplexLU ((Record (: X (Vector (Float))) (: objVal (Float)) (: Basis (List (Integer)))) (Matrix (Float)) (Vector (Float)) (Vector (Float)) (List (Integer))))): Not documented!!!!
--->-->LinearProgramming((luDecompose ((Record (: L (Matrix (Float))) (: U (Matrix (Float))) (: P (Matrix (Float)))) (Matrix (Float))))): Not documented!!!!
--->-->LinearProgramming((luSolve ((Record (: X (Vector (Float))) (: P (Matrix (Float)))) (Matrix (Float)) (Vector (Float))))): Not documented!!!!
; compiling file "/var/aw/var/LatexWiki/LINPROG.NRLIB/LINPROG.lsp" (written 16 JUL 2018 07:35:26 PM):
; /var/aw/var/LatexWiki/LINPROG.NRLIB/LINPROG.fasl written
; compilation finished in 0:00:00.285
------------------------------------------------------------------------
LinearProgramming is now explicitly exposed in frame initial
LinearProgramming will be automatically loaded when needed from
/var/aw/var/LatexWiki/LINPROG.NRLIB/LINPROG

Test flavours

fricas
-- Unittest .....: Package LinearProgramming
Author .......: Kurt Pagani -- Date .........: Sun Jun 17 02:14:10 CEST 2018
fricas
)set break resume

fricas
)expose UnittestCount UnittestAux Unittest
UnittestCount is now explicitly exposed in frame initial
UnittestAux is now explicitly exposed in frame initial
Unittest is now explicitly exposed in frame initial
--)co linprog
fricas
)expose LINPROG
LinearProgramming is already explicitly exposed in frame initial
-- Test Suite
testsuite "Package: LinearProgramming"
All user variables and function definitions have been cleared.
WARNING: string for testsuite should have less than 15 characters!
Type: Void
fricas
-- Cases
testcase "revisedSimplex"
All user variables and function definitions have been cleared.
Type: Void
fricas
-- testEquals("", "")
A:=matrix [[0.7,1.0,1.0,0.0,0.0,0.0],[0.5,(5/6)::Float,0.0,1.0,0.0,0.0],
[1.0,(2/3)::Float,0.0,0.0,1.0,0.0],[0.1,0.25,0.0,0.0,0.0,1.0]]
 (1)
Type: Matrix(Float)
fricas
c:=vector [-10.0,-9.0,0.0,0.0,0.0,0.0]
 (2)
Type: Vector(Float)
fricas
b:=vector [630.0,600.0,708.0,135.0]
 (3)
Type: Vector(Float)
fricas
B:=[3,4,5,6]
 (4)
Type: List(PositiveInteger?)
fricas
X1:Vector Float:= vector [540.0,252.0,0.0,120.0,0.0,18.0]
 (5)
Type: Vector(Float)
fricas
objVal1:Float:= - 7668.0
 (6)
Type: Float
fricas
Basis1:List Integer:= [2,4,1,6]
 (7)
Type: List(Integer)
fricas
eps:=0.0000000001
 (8)
Type: Float
fricas
floatNull(v:Vector Float):Boolean == sqrt dot(v,v) <= eps
Function declaration floatNull : Vector(Float) -> Boolean has been
added to workspace.
Type: Void
fricas
r1:=revisedSimplex( A, b, c, B)
 (9)
Type: Record(X: Vector(Float),objVal: Float,Basis: List(Integer))
fricas
testEquals("floatNull(r1.X-X1)", "true")
fricas
Compiling function floatNull with type Vector(Float) -> Boolean
Type: Void
fricas
testEquals("abs(r1.objVal-objVal1)<eps", "true")
Type: Void
fricas
testEquals("r1.Basis", "Basis1")
Type: Void
fricas
r1lu:=revisedSimplexLU( A, b, c, B)
 (10)
Type: Record(X: Vector(Float),objVal: Float,Basis: List(Integer))
fricas
testEquals("floatNull(r1lu.X-X1)", "true")
Type: Void
fricas
testEquals("abs(r1lu.objVal-objVal1)<eps", "true")
Type: Void
fricas
testEquals("r1lu.Basis", "Basis1")
Type: Void
fricas
--
---
A2 := matrix [[0,1,-1,0,0,0], [1,1,0,-1,0,0],
[-0.5,1,0,0,1,0],[-1,1,0,0,0,1]]
 (11)
Type: Matrix(Float)
fricas
c2 := vector [2,-3,0,0,0,0]
 (12)
Type: Vector(Integer)
fricas
b2 := vector [1,2,8,6]
 (13)
Type: Vector(PositiveInteger?)
fricas
B2 := [1,2,5,6]
 (14)
Type: List(PositiveInteger?)
fricas
X2:= [4.0,10.0,9.0,12.0,0.0,0.0]
 (15)
Type: List(Float)
fricas
objVal2:= - 22.0
 (16)
Type: Float
fricas
Basis2:= [3,2,1,4]
 (17)
Type: List(PositiveInteger?)
fricas
r2:=revisedSimplex( A2, b2, c2, B2)
 (18)
Type: Record(X: Vector(Float),objVal: Float,Basis: List(Integer))
fricas
testEquals("floatNull(r2.X-X2)", "true")
Type: Void
fricas
testEquals("abs(r2.objVal-objVal2)<eps", "true")
Type: Void
fricas
testEquals("r2.Basis", "Basis2")
Type: Void
fricas
r2lu:=revisedSimplexLU( A2, b2, c2, B2)
 (19)
Type: Record(X: Vector(Float),objVal: Float,Basis: List(Integer))
fricas
testEquals("floatNull(r2lu.X-X2)", "true")
Type: Void
fricas
testEquals("abs(r2lu.objVal-objVal2)<eps", "true")
Type: Void
fricas
testEquals("r2lu.Basis", "Basis2")
Type: Void
fricas
--
A3 := matrix [[0.4,0.5,1,0,0],[0,0.2,0,1,0],[0.6,0.3,0,0,1]]
 (20)
Type: Matrix(Float)
fricas
c3 := vector [-40,-30,0,0,0]
 (21)
Type: Vector(Integer)
fricas
b3 := vector [20,5,21]
 (22)
Type: Vector(PositiveInteger?)
fricas
B3 := [3,4,5]
 (23)
Type: List(PositiveInteger?)
fricas
X3:= [25.0,20.0,0.0,1.0,0.0]
 (24)
Type: List(Float)
fricas
objVal3:= - 1600.0
 (25)
Type: Float
fricas
Basis3:= [2,4,1]
 (26)
Type: List(PositiveInteger?)
fricas
r3:=revisedSimplex(A3, b3, c3, B3)
 (27)
Type: Record(X: Vector(Float),objVal: Float,Basis: List(Integer))
fricas
testEquals("floatNull(r3.X-X3)", "true")
Type: Void
fricas
testEquals("abs(r3.objVal-objVal3)<eps", "true")
Type: Void
fricas
testEquals("r3.Basis", "Basis3")
Type: Void
fricas
r3lu:=revisedSimplexLU(A3, b3, c3, B3)
 (28)
Type: Record(X: Vector(Float),objVal: Float,Basis: List(Integer))
fricas
testEquals("floatNull(r3lu.X-X3)", "true")
Type: Void
fricas
testEquals("abs(r3lu.objVal-objVal3)<eps", "true")
Type: Void
fricas
testEquals("r3lu.Basis", "Basis3")
Type: Void
fricas
--
A4 := matrix [[1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],_
[0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],_
[0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],_
[0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0],_
[3, 0, 6, 0, 5, 0, 7, 0, 1, 0, 0, 0, 0, 0, 1, 0],_
[0, 2, 0, 4, 0, 10, 0, 4, 0, 1, 0, 0, 0, 0, 0, 1]]
 (29)
Type: Matrix(NonNegativeInteger?)
fricas
b4 := vector [1, 1, 1, 1, 13, 10]
 (30)
Type: Vector(PositiveInteger?)
fricas
c4 := vector [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]
 (31)
Type: Vector(NonNegativeInteger?)
fricas
B4 := [11, 12, 13, 14, 15, 16]
 (32)
Type: List(PositiveInteger?)
fricas
X4:= [1.0, 0.0, 0.25, 0.7500000000_0000000001, 0.3, 0.7, 1.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 (33)
Type: List(Float)
fricas
objVal4:= 0.0
 (34)
Type: Float
fricas
Basis4:= [1,3,6,5,7,4]
 (35)
Type: List(PositiveInteger?)
fricas
r4 := revisedSimplex( A4, b4, c4, B4)
 (36)
Type: Record(X: Vector(Float),objVal: Float,Basis: List(Integer))
fricas
testEquals("floatNull(r4.X-X4)", "true")
Type: Void
fricas
testEquals("abs(r4.objVal-objVal4)<eps", "true")
Type: Void
fricas
testEquals("r4.Basis", "Basis4")
Type: Void
fricas
r4lu := revisedSimplexLU( A4, b4, c4, B4)
 (37)
Type: Record(X: Vector(Float),objVal: Float,Basis: List(Integer))
fricas
testEquals("floatNull(r4lu.X-X4)", "true")
Type: Void
fricas
testEquals("abs(r4lu.objVal-objVal4)<eps", "true")
Type: Void
fricas
testEquals("r4lu.Basis", "Basis4")
Type: Void
fricas
--
A5 := matrix [[1, -1, 1, 0, 0], [4, 9, 0, 1, 0], [-2, 4, 0, 0, 1]]
 (38)
Type: Matrix(Integer)
fricas
c5 := vector [-1, -1, 0, 0, 0]
 (39)
Type: Vector(Integer)
fricas
b5 := vector [2, 18, 4]
 (40)
Type: Vector(PositiveInteger?)
fricas
B5 := vector [3, 4, 5]
 (41)
Type: Vector(PositiveInteger?)
fricas
X5 := [2.7692307692_307692308, 0.7692307692_3076923077, 0.0, 0.0,
6.4615384615_384615385]
 (42)
Type: List(Float)
fricas
objVal5:= - 3.5384615384_615384615
 (43)
Type: Float
fricas
Basis5:= [1,2,5]
 (44)
Type: List(PositiveInteger?)
fricas
r5 := revisedSimplex( A5, b5, c5, B5)
 (45)
Type: Record(X: Vector(Float),objVal: Float,Basis: List(Integer))
fricas
testEquals("floatNull(r5.X-X5)", "true")
Type: Void
fricas
testEquals("abs(r5.objVal-objVal5)<eps", "true")
Type: Void
fricas
testEquals("r5.Basis", "Basis5")
Type: Void
fricas
r5lu := revisedSimplexLU( A5, b5, c5, B5)
 (46)
Type: Record(X: Vector(Float),objVal: Float,Basis: List(Integer))
fricas
testEquals("floatNull(r5lu.X-X5)", "true")
Type: Void
fricas
testEquals("abs(r5lu.objVal-objVal5)<eps", "true")
Type: Void
fricas
testEquals("r5lu.Basis", "Basis5")
Type: Void
fricas
-- Results/Statistics
fricas
)set output algebra on

fricas
)version
Value = "FriCAS 1.3.4 compiled at Wed Jun 27 14:53:01 UTC 2018"
fricas
)lisp (lisp-implementation-type)
Value = "SBCL"
fricas
)lisp (lisp-implementation-version)
Value = "1.1.1"
statistics()
=============================================================================
General WARNINGS:
* do not use ')clear completely' before having used 'statistics()'
It clears the statistics without warning!
* do not forget to pass the arguments of the testXxxx functions as Strings!
Otherwise, the test will fail and statistics() will not notice!
* testLibraryError does not prevent FriCAS from aborting the current block.
Thus, if a block contains other test functions, they will not be executed
and statistics() will not notice!
=============================================================================
Testsuite: Package: LinearProgramming
failed (total): 0 (1)
=============================================================================
testsuite | testcases: failed (total) | tests: failed (total)
Package: LinearProgramming    0     (1)               0    (30)
=============================================================================
File summary.
unexpected failures: 0
expected failures: 0
unexpected passes: 0
total tests: 30
Type: Void

 Subject:   Be Bold !! ( 14 subscribers )