(defmacro vector-push (n a fp) `(cond ((>= ,fp (array-dimension ,a 0)) nil) (t (setf (aref ,a ,fp) ,n) (incf ,fp) (- ,fp 1)))) (proclaim '(optimize (speed 3) (safety 0) (space 0) (compilation-speed 0))) ;; Simplex Algorithm and its demonstration ;; Bruno Haible 25.09.1986, 30.5.1991, 22.4.1992 ;; Common Lisp version ;; English version ;; Copyright (c) Bruno Haible 1991, 1992.. ;; This file may be copied under the terms of the GNU General Public License. (provide 'simplex) ; type of a number: here rational numbers (deftype Zahl () 'float) ; abbreviation: R ; Should this be changed to non-exact numbers (floats), uncomment the ; portion #| ... |# of code and eliminate the calls to the function rational. (deftype Zahl-Vektor () '(simple-array Zahl 1)) ; abbreviation: R^n (deftype Zahl-Matrix () '(simple-array Zahl 2)) ; abbreviation: R^mxn ; Notation: ; indices 1-based in mathematical theory, 0-based in this implementation ; A` the transpose of A ; x[i] (aref x i) ; A[i,j] (aref A i j) ; u`*v = sum(i>=1, u[i]*v[i]) ; Input: A, an mxn matrix of numbers, ; b, an m vector of numbers, ; c, an n vector of numbers. ; Solves a canonical optimisation problem: ; ( C ) Among the y in R^n with A*y=b, y>=0 find those with = Min ; ( C* ) Among the z in R^m with x=A`*z+c>=0 find those with = Min ; See: Koulikov, Alg`ebre et th'eorie des nombres, Chap. 6.1. ; Output: 1st value: flag if solvable, NIL or T, ; if solvable: ; 2nd value: solution of ( C ), the linear problem LP: ; list containing ; - the solution y ; - the optimal value of ; 3rd value: solution of ( C* ), the dual problem DP: ; list containing ; - the solution x ; - the constant terms Z for the z's ; - the dependencies ZZ between the z's ; - the optimal value of ; The general solution has the following form: ; If ZZ[i,i]=T, the variable z_i is unrestricted. ; Otherwise ZZ[i,i]=NIL, and the variable z_i is calculated by ; z_i = Z[i] + sum(j=1,...,m with Z[j,j]=T , ZZ[i,j]*z_j) . ; 4th value: flag if (T) the given solutions make up the whole ; set of solutions or (NIL) there may be another solution. ; In either case the set of solutions is convex. ; Method: ; Method of Dantzig, see Koulikov book; ; enhancement for degenerate case following Bronstein/Semendjajew. ; [Thanks to all these russian mathematicians!] (defun simplex (A b c) ; check input: (let ((m (length b)) (n (length c))) (declare (type fixnum m n)) (assert (typep A `(array * (,m ,n)))) ; m = height of the matrix we are working in. ; n = width of the matrix we are working in. ; initialisation: (let ((AA (let ((h (make-array `(,m ,n)))) (dotimes (i m) (dotimes (j n) (setf (aref h i j) (aref A i j)) ) ) h ) ) (AB (let ((h (make-array m))) (dotimes (i m) (setf (aref h i) (elt b i))) h ) ) (AC (let ((h (make-array n))) (dotimes (j n) (setf (aref h j) (elt c j))) h ) ) (AD 0) ; AA in R^mxn : the tableau we work in ; AB in R^m : the right margin ; AC in R^n : the bottom margin ; AD in R : the bottom right corner (Extend nil) (AAE nil) ; Extend in Boolean : flag if the additional columns (degenerate case) are being used ; AAE in R^mxm : additional columns for the degenerate case (Zeile (let ((h (make-array m))) (dotimes (i m) (setf (aref h i) (cons 'ZO i))) h ) ) (Spalte (let ((h (make-array n))) (dotimes (j n) (setf (aref h j) (cons 'XY j))) h ) ) ; Zeile in Cons^m : left labels ; Spalte in Cons^n : top labels ; ; The tableau: ; ; | Zeile[0] Zeile[1] ... Zeile[n-1] | ; ------------+-------------------------------------+--------- ; Spalte[0] | AA[0,0] AA[0,1] ... AA[0,n-1] | AB[0] ; Spalte[1] | AA[1,0] AA[1,1] ... AA[1,n-1] | AB[1] ; ... | ... ... ... | ... ; Spalte[m-1] | AA[m-1,0] AA[m-1,1] ... AA[m-1,n-1] | AB[m-1] ; ------------+-------------------------------------+--------- ; | AC[0] AC[1] ... AC[n-1] | AD ; ; Labeling of columns: 0 or Y ; ... ... ; Z X ; ; Labeling of rows: Z ... 0 or X ... -Y. ; ; for all i=0,...,m-1: ; sum(j=0,...,n-1; AA[i,j]*(0 or Y[Spalte[j].Nr])) - AB[i] = ; = (0 or -Y[Zeile[i].Nr]) ; for all j=0,...,n-1: ; sum(i=0,...,m-1; AA[i,j]*(Z[Zeile[i].Nr] or X[Zeile[i].Nr])) + AC[j] = ; = (Z[Spalte[j].Nr] or X[Spalte[j].Nr]) ; sum(j=1,...,N; AC[j]*(0 or Y[Spalte[j].Nr])) - AD = ; sum(i=1,...,M; AB[i]*(Z[Zeile[i].Nr] or X[Zeile[i].Nr])) + AD = ; These are to be considered as equations in the unknowns X,Y,Z. ; ; The additional columns - if present - are added at the right. ; (ZZ (let ((h (make-array `(,m ,m) :initial-element 0))) (dotimes (i m) (setf (aref h i i) nil)) h ) ) (Z (make-array m)) (Y (make-array n)) (X (make-array n)) ) (declare (type Zahl-Matrix AA) (type Zahl-Vektor AB AC) (type Zahl AD)) (flet ((pivot (k l) (declare (type fixnum k l)) ; pivots the tableau around the element AA[k,l] with 0<=k h hmax) (setq hmax h l j)) ) ) ) (if l ; AA[i,l] was the maximal element w.r.t. absolute value (progn (vector-push i not-elbar not-elbar-fp) (pivot i l) ) ; trivial row (if (zerop (aref AB i)) ; Keep this dummy line; it will not change since we will ; pivot only around XY columns and these columns have a 0 ; in row i. ; Nonzero elements are only in ZO columns, and these ; will not change. (vector-push i elbar elbar-fp) ; this row makes LP inconsistent -> unsolvable (return-from simplex NIL) ) ) ) ) ; Mark the eliminatable row in the ZZ matrix: ; The non-eliminatable rows have swapped positions with the X ; such that we have Spalte[(cdr Zeile[i])] = (ZO . i) if row i ; is non-eliminatable ( <==> (car Zeile[i]) = XY ). (dotimes (i0h elbar-fp) (let ((i0 (aref elbar i0h))) (setf (aref Z i0) 0) ; we must set Z[i0]=0 to make the other Z's correct! (setf (aref ZZ i0 i0) T) ; z_i0 is unrestricted (dotimes (ih not-elbar-fp) (let* ((i (aref not-elbar ih)) ; we must have (car Zeile[i])=XY (j (cdr (aref Zeile i)))) ; column with which row i was pivoted (setf (aref ZZ i i0) (aref AA i0 j)) ) ) ) ) ; delete rows: (uses that every not-elbar[i]>=i) (dotimes (ih not-elbar-fp) (let ((i (aref not-elbar ih))) (unless (eql ih i) (dotimes (j n) (setf (aref AA ih j) (aref AA i j))) (setf (aref AB ih) (aref AB i)) ) ) ) (setq m not-elbar-fp) ; new number of rows = number of the ZO columns ) ; sort columns: bring XY to the left, ZO to the right. ; This is used at the end to calculate the non-eliminatable z from the X. (let ((l 0) ; left column (r (1- n))) ; right column (loop (unless (< l r) (return)) (cond ((eq (car (aref Spalte l)) 'XY) (incf l)) ((eq (car (aref Spalte r)) 'ZO) (decf r)) (t ; swap columns r and l (dotimes (i m) (rotatef (aref AA i l) (aref AA i r))) (rotatef (aref AC l) (aref AC r)) (rotatef (aref Spalte l) (aref Spalte r)) ) ) ) ) ; hide these M columns from pivoting: (setq n (- n m)) ; The elements AA[0..m-1,n..n+m-1], AC[n..n+m-1], Spalte[n,n+m-1] ; will only be used again at the end of phase 6. (let ((Zeile_save (copy-seq Zeile))) (flet ((SuchePivotZeile (AWZM l AWZM-fp) ; For a choice set AWZM of rows and a column l choose the ; row k such that: ; We assume that ; for i in AWZM we have 0<=i=0 und AA[i,l]>0. ; Among the i in AWZM, k is the one for which the quotient ; AB[i]/AA[i,l] is minimal. Should this quotient be =0, or ; if Extend=true, then afterwards Extend=true, and the vector ; (AB[i]/AA[i,l], AAE[i,1]/AA[i,l], ..., AAE[i,m]/AA[i,l]) ; has been minimized (lexicographically) among all i in AWZM. ; If AWZM is empty, NIL is returned. (if (eql AWZM-fp 0) NIL (let (k) (unless Extend ; try to choose k (let (hmax) (dotimes (ih AWZM-fp) (let* ((i (aref AWZM ih)) (h (/ (aref AB i) (aref AA i l)))) (when (or (eql ih 0) (< h hmax)) (setq hmax h k i)) ) ) (when (zerop hmax) ; degenerate case (setq Extend T) (setq AAE (make-array `(,m ,m) :initial-element 0)) (dotimes (i m) (setf (aref AAE i i) 1)) ) ) ) (when Extend ; The degenerate case has already been active or has ; emerged now (then the old k may be bad). (let (hmax hmaxe) (dotimes (ih AWZM-fp) (let ((i (aref AWZM ih))) (let ((h (/ (aref AA i l))) (he (make-array m))) (dotimes (j m) (setf (aref he j) (* (aref AAE i j) h))) (setq h (* (aref AB i) h)) ; (h,he[1],...,he[M]) is the vector of quotients (when (or (eql ih 0) ; instead of (< h hmax) now lexicographic comparison ; (h,he[1],...,he[M]) < (hmax,hmaxe[1],...,hmaxe[M]) : (or (< h hmax) (and (= h hmax) (dotimes (ie m NIL) (let* ((he_ie (aref he ie)) (hmaxe_ie (aref hmaxe ie))) (when (< he_ie hmaxe_ie) (return T)) (when (> he_ie hmaxe_ie) (return NIL)) ) ) ) ) ) (setq hmax h hmaxe he k i) ) ) ) ) ) ) k )) ) ) (let ((PZM (make-array m :fill-pointer 0)) (PZM-fp 0) (p nil)) ; PZM = set of the i with AB[i]>=0 ; p = the last i for which AB[i] was being maximized (loop ; fill column AB with positive numbers ; throw away roundoff errors: (dotimes (ih PZM-fp) (let ((i (aref PZM ih))) (when (minusp (aref AB i)) (setf (aref AB i) 0)) ) ) (let ((NZM (make-array m :fill-pointer 0)) (NZM-fp 0)) ; NZM = set of the i with AB[i]<0 ; recalculate PZM and NZM: (let ((old-PZM-count PZM-fp)) ; old cardinality of PZM (setf PZM-fp 0) (dotimes (i m) (if (>= (aref AB i) 0) (vector-push i PZM PZM-fp) (vector-push i NZM NZM-fp) ) ) ; delete the additional columns if PZM really grew ; and the degeneracy perhaps disappeared: (when (> PZM-fp old-PZM-count) (setq Extend nil) (setq p nil) ) ; otherwise PZM remained unchanged, and AB[p]<0 still holds. ) (when (eql NZM-fp 0) (return)) ; every AB[i]>=0 ? (if p ; use the last p. (when (dotimes (j n t) (when (< (aref AA p j) 0) (return nil))) ; every AA[p,j]>=0 but AB[p]<0 ; ==> row makes LP inconsistent ==> unsolvable (return-from simplex NIL) ) ; choose p: p := the i among those with AB[i]<0 ; for which the number of AA[i,j]>=0 is maximal. (let ((countmax -1)) (dotimes (ih NZM-fp) (let ((i (aref NZM ih)) (count 0)) (dotimes (j n) (when (>= (aref AA i j) 0) (incf count))) (when (> count countmax) (setq countmax count p i)) ) ) (when (eql countmax n) ; every AA[p,j]>=0 but AB[p]<0 ; ==> row makes LP inconsistent ==> unsolvable (return-from simplex NIL) ) ) ) ; Now AB[p]<0, and there is a j with AA[p,j]<0. ; Choose l: maximal abs(AA[p,j]) among the j with AA[p,j]<0. (let ((hmin 0) (l nil)) (dotimes (j n) (let ((h (aref AA p j))) (when (< h hmin) (setq hmin h l j)) ) ) ; build AWZM: (let ((AWZM (make-array m :fill-pointer 0)) (AWZM-fp 0)) (dotimes (ih PZM-fp) (let ((i (aref PZM ih))) (when (> (aref AA i l) 0) (vector-push i AWZM AWZM-fp)) ) ) (let ((k (SuchePivotZeile AWZM l AWZM-fp))) (if (null k) ; Pivoting around AA[p,l] lets PZM grow at least ; by the element p because: ; for i in PZM we have AB[i]>=0 and ; (because AZWM={}) also AA[i,l]<=0, ; and then afterwards ; AB[i]=AB[i]-AB[p]*AA[i,l]/AA[p,l] >= 0, ; >=0 <0 <=0 <0 ; that is i in PZM again. ; But afterwards we also have ; AB[p]=AB[p]/AA[p,l] (<0/<0) >0, ; therefore p in PZM. (progn (setq Extend nil) ; additional columns aren't needed any more (pivot p l) ) ; pivot row k chosen, ready for pivoting. ; Pivoting around AA[k,l] does not shrink PZM ; because: ; for i in PZM we have AB[i]>=0. ; If AA[i,l]<=0, we have afterwards ; AB[i]=AB[i]-AB[k]*AA[i,l]/AA[k,l] >=0. ; >=0 >=0 <=0 >0 ; But if AA[i,l]>0, then for every i/=k afterwards ; AB[i]=AA[i,l]*(AB[i]/AA[i,l]-AB[k]/AA[k,l]) >=0 ; >0 >=0 because of the choice of k and i in AWZM ; and for i=k afterwards AB[k]=AB[k]/AA[k,l] (>=0/>0) >=0. ; Therefore afterwards always AB[i]>=0, i.e. i in PZM. (pivot k l) ) ) ) ) ) ) ) ; Now Extend=false, since (fill-pointer PZM) must just have ; grown, reaching m. ; From now on every AB[i]>=0, and this property remains. (loop ; search an l with AC[l]<0 : (let (l) (when (dotimes (j n T) (when (< (aref AC j) 0) (setq l j) (return NIL)) ) (return) ; every AC[j]>=0 ==> solvable ) ; AWZM := set of the i with AA[i,l]>0 : (let ((AWZM (make-array m :fill-pointer 0)) (AWZM-fp 0)) (dotimes (i m) (when (> (aref AA i l) 0) (vector-push i AWZM AWZM-fp)) ) (let ((k (SuchePivotZeile AWZM l AWZM-fp))) (if (null k) ; every AA[i,l]<=0 and AC[l]<0 ==> column makes DP inconsistent (return-from simplex NIL) (pivot k l) ; still every AB[i]>=0. ; AD:=AD-AB[k]*AC[l]/AA[k,l] >= AD is not lowered. ; >=0 <0 >0 ) ) ) ) ) ; Solvable! Build solution: (let ((complete t)) (dotimes (i m) (let ((s (aref AB i)) (index (cdr (aref Zeile i)))) (setf (aref X index) 0 (aref Y index) s) (setq complete (and complete (> s 0))) ) ) (dotimes (j n) (let ((s (aref AC j)) (index (cdr (aref Spalte j)))) (setf (aref X index) s (aref Y index) 0) (setq complete (and complete (> s 0))) ) ) ; The non-eliminatable z values are calculated from the hidden ; parts of AA and AC and the values of X and Y : (do ((j n (1+ j))) ((>= j (+ n m))) (let ((s (aref AC j))) (dotimes (i m) (setq s (+ s (* (aref AA i j) (aref X (cdr (aref Zeile_save i)))))) ) (setf (aref Z (cdr (aref Spalte j))) s) ) ) (values T (list Y (- AD)) (list X Z ZZ AD) complete) ) ) ) ) ) ) )