Simplex Algorithm
11-21-2023, 11:20 PM (This post was last modified: 11-23-2023 07:48 PM by Albert Chan.)
Post: #41
 Albert Chan Senior Member Posts: 2,611 Joined: Jul 2018
RE: Simplex Algorithm
(11-21-2023 05:58 PM)Albert Chan Wrote:  With vars moved to the end, user see [min(c) ... vars]

I tracked number of pivot1 calls inside simplex_core(), simplex algorithm is very efficient.

simplex test_cases pivot1 calls.

1: 3+3 = 6      // use simplex2 wrapper
2: 3+3 = 6
3: 3+2 = 5
4: 2+3 = 5
5: 3+3 = 6
6: 2+2 = 4
7: 2+2 = 4      // simplex_int use simplex2 once
8: 2+3 = 5
9: 2+11 = 13
10: 1+1 = 2
11: 2+10 = 12

1st number is simplex2 len(bv) pivot1 calls, to make all pivots=1.
2nd number is simplex_core() total pivot1 calls, to reach solution.

Here is version 5, added all posted ideas since 4
Code:
#cas simplex(a):= BEGIN LOCAL c,v,colJ,I,J,bv,n,(art:=0); IF dim(dim(a))!=2 THEN a,art:=a END n:=art; c,v:=dim(a).-1;  /* constraints, variables */ bv:=makelist(c); /* pivots */ FOR J FROM v DOWNTO 1 DO  IF n>0 AND J<=v-art THEN error(1,a) END  colJ:=col(a,J)[1..c];  IF (I:=contains(colJ,1))==0 THEN continue END  IF bv[I] or l1norm(colJ[I]:=0) THEN continue END  bv[I]:=J;  IF c > art-(n-=1) THEN continue END  RETURN simplex_core(addartrow(a,art),bv,v,art,0) END; error(2,a)  /* too few pivots */ END; simplex2(a,bv,v,art,ign):= BEGIN  LOCAL I;  FOR I FROM 1 to len(bv) DO a:=pivot1(a,I,bv[I]) END  RETURN simplex_core(a,bv,v,art,ign) END; simplex2(a,bv,v,art,ign):= BEGIN  LOCAL I;  FOR I FROM 1 to len(bv) DO a:=pivot1(a,I,bv[I]) END  RETURN simplex_core(a,bv,v,art,ign) END; simplex_core(a,bv,v,art,ign):= BEGIN LOCAL a2,b,c,d,bmin,cmin,n,l1,l2,I,J,P:=0; d:=dim(a); c:=SUPPRESS(row(a,d(1)),d(2)-ign,d(2)); cmin:=MIN(c); n:=1+(art>0); b:=SUPPRESS(col(a,0),d(1)-n+1,d(1)); bmin:=MIN(b); // Step 2: If all ci,bi≥0 go to Step 6 WHILE cmin<0 OR bmin<0 DO // Step 3: For any ci < 0 check for unbounded solution  FOR I FROM 1 TO (d(2)-1-ign) DO   IF c(I) < 0 AND MAX(col(a,I))<=0 THEN // UNBOUNDED SOLUTION    RETURN [-inf,a,bv,P,[]];   END;  END; // Step 4: Determine bv to enter and exit  J:=POS(c,cmin); // xin  I:=minRatio(col(a,J),col(a,0),n); // xout // Step 5: Canonical reduction with new bv  IF bv[I]!=J THEN   bv[I]:=J;   a:=pivot1(a,I,J); P+=1;   c:=SUPPRESS(row(a,d(1)),d(2)-ign,d(2));   cmin:=MIN(c);  END;  IF cmin>=0 THEN   IF art>0 THEN    IF a(d(1),d(2))!=0 THEN // w≠0? NO SOLUTION / EMPTY FEASIBLE REGION     RETURN ["Bad Row "+d(1),a,bv,P,[]];    END;    a:=delrows(a,d[1]);    d[1]-=(n:=1);    ign:=art;    art:=0;    c:=SUPPRESS(row(a,d(1)),d(2)-ign,d(2));    cmin:=MIN(c);   ELSE    b:=SUPPRESS(col(a,0),d(1)-n+1,d(1));    bmin:=MIN(b); // IF ∃bi<0, use Dual Simplex    IF bmin<0 THEN     FOR I FROM 1 TO len(b) DO // For any bi < 0 check for unbounded/no solution      IF b(I)<0 AND MIN(SUPPRESS(row(a,I),d(2)-ign,d(2)))>=0 THEN       RETURN ["Bad Row "+I,a,bv,P,[]];      END;     END;     I:=POS(b,bmin);     J:=minRatio(-a[I],c,0); /* max ratio */     bv[I]:=J;     a:=pivot1(a,I,J); P+=1;     c:=SUPPRESS(row(a,d(1)),d(2)-ign,d(2));     cmin:=MIN(c);    END;   END;  END; END; // Step 6: ci,bi≥0 optimum achieved l1:=sol(a,v); a2:=a; c:=SUPPRESS(row(a2,d(1)),d(2)-ign,d(2)); FOR J FROM 1 TO len(c) DO  IF c(J) or contains(bv,J) THEN continue END  IF (b:=abs(trn(col(a2,J))))==1 or b==0 THEN continue END  I:=minRatio(col(a2,J),col(a2,0),1);  a2:=pivot1(a2,I,J); P+=1;  l2:=sol(a2,v); // IFINITE SOLUTIONS  IF l1==l2 THEN BREAK END  RETURN [-a[0][0],a,bv,P,transpose(l1,l2)]; END; // UNIQUE SOLUTION RETURN [-a[0][0],a,bv,P,l1]; END; sol(a,v):= BEGIN  LOCAL x,b,c, j,k;  x:=makelist(v); /* all zeroes */  b:=col(a,0);    /* last column */  FOR j FROM 1 TO v DO   IF abs(trn(c:=col(a,j)))!=1 THEN continue END   IF (k:=contains(c,1))==0 THEN continue END   IF b[k] THEN x[j]:=b[k]; b[k]:=0 END  END;  RETURN x; END; minRatio(a,b,n):=  BEGIN  LOCAL i,ratio, (j:=0),(minratio:=inf);  FOR i FROM 1 TO len(b)-n DO   IF a[i]<=0 THEN continue END   ratio:=b[i]/a[i];   IF minratio>ratio THEN minratio:=ratio; j:=i END  END;  RETURN j END; // first attempt of Gomory's plane cutting algorithm simplex_int(a,bv,v,art,ign,integers):= BEGIN  LOCAL b,d,x,n; // Step 1: ignore integer restriction, apply simplex  a:=simplex2(a,bv,v,art,ign);  ign:=art;  WHILE 1 DO   IF len(x:=a[0])==0 THEN return a END // no solution   IF len(x[1])>1 THEN return a END     // inf solution      FOR I FROM 1 to len(integers) DO    IF type(x[integers[I]])!=DOM_INT THEN break END   END;   IF I>len(integers) THEN return a END // all integers // Step 2: Detetmine constraint with largest fractional part   a,bv,n := a[2..4];   d:=dim(a);   b:=col(a,d[2]);   x := b[1 .. d[1]-1];   x -= floor(x);   x := a[POS(x, max(x))]; // Create new constraint   x:=floor(x)-x;   I:=len(x)-ign;   x:=extend(x[1..I-1], [1], x[I..0]);   bv:=append(bv,d(2)-ign); /* pivot on 1 */    // Step 3: solve with new constraint   REDIM(a,{d(1),d(2)-1});   ADDCOL(a,b-b,d(2)-ign);   ADDCOL(a,b,d(2)+1-ign);   ADDROW(a,x,d(1));   a:=simplex_core(a,bv,v,0,ign);   a[4]+=n;  // pivot1 calls  END; END; simplex_le(a) := BEGIN  LOCAL b, p, v,(p:=[]),(d:=inf);  IF dim(dim(a))!=2 THEN a,d:=a END  IF abs(d)<len(a) THEN a,p:=remove_eq_c(a,d) END  b := dim(a);  a := transpose(extend(a, identity(b[1])));  a[0] := a[b[2]];  a := transpose(delrows(a, b[2] .. b[2]+len(p)));  IF d==0 THEN return a END  v := sum(b)-2 - len(p);         /* variables */  p := extend(p,range(b[2],v+1)); /* pivots */  IF d>0 THEN a[0] := -a[0] END   /* maximize */  a := simplex_core(a,p,v,0,0);  IF d>0 THEN a[1] := -a[1] END   /* undo flip */  return a; END; remove_eq_c(a, n) := BEGIN  LOCAL I,J,c := len(a[1]);  LOCAL p := makelist(n:=abs(n)); /* max pivots */  FOR I FROM n DOWNTO 1 DO   p[I] := J := pivot_index(suppress(a[I],c));   IF J==0 and a[I][c] THEN error(3,a) END   a := J ? pivot1(a,I,J) : delrows(a,I);  END  return a, remove(0,p); END; pivot_index(a) := BEGIN  LOCAL b:=contains((a:=abs(a)),1);  IF b THEN RETURN b END  IF (b:=max(a)) THEN b:=contains(a,b) END  RETURN b; END; pivot1(m,I,J):= pivot((m[I]:=m[I]/m[I][J]),I,J); artrow(col,art):=append(makelist(k->k>col-art,2,col),0); addartrow(a,art):=  /* art vars are valid pivots */ BEGIN  local c:=len(a[1]);  IF art<=0 THEN RETURN a END  append(a,artrow(c,art)-sum(a[POS(col(a,J),1)],J,c-art,c-1))   END; #end

simplex_core() is fast, but assume all supplied pivots are good (unit column).
It now returns [min(c), matrix, bv, pivot1_calls, vars]
On HP Prime, most likely, is showed as [min(c) ... vars]

simplex2() is the safer simplex_core(), go thru all the pivots, to make sure all 1's.

simplex() is also very safe, it confirmed a unit column before using it as pivot.
And because it only use true pivot, it does not need pivot1() calls.

It is not simply adding artrow to a, but also to make art variables available to use as pivots.
Note: if art <= 0 it just return a. i.e. no art row is added.

Example:
> m:=[[1,2,0,1,1,0,0,20],[2,1,1,0,0,1,0,10],[-1,4,-2,3,0,0,1,40],[1,4,3,2,0,0,0,0]]

[35 ... [5,0,0,15,0,0,0]]

This is equivalent to simplex2 with plain artrow added, but faster, without pivot1() calls.

> simplex2(append(m,artrow(8,3)),[5,6,7],7,3,0)

[35 ... [5,0,0,15,0,0,0]]

My previous version hide this detail, with artrow added inside simplex_core.
But, to make core restartable at any stage, we need to expose this artrow.