-- TableauImpl.mesa
-- Last Edited by: Gnelson, November 21, 1983 4:59 pm
-- Last Edited by: Beach, January 12, 1984 1:35 pm

DIRECTORY Tableau, Real;

TableauImpl: PROGRAM IMPORTS Tableau EXPORTS Tableau = {
OPEN Tableau;

matrix: PUBLIC REF ARRAY[1 .. nrows] OF ARRAY[1 .. ncols] OF REAL _
  NEW[ARRAY[1 .. nrows] OF ARRAY[1 .. ncols] OF REAL];

n, m: PUBLIC INT;

dcol: PUBLIC INT;

unit: PUBLIC Variable _ NewVariable[];

Init: PUBLIC PROC = {n _ 0; m _ 1; unit.i _ 0; unit.j _ 1; satisfiable _ TRUE};

satisfiable: PUBLIC BOOL;

NewVariable: PUBLIC PROC RETURNS [Variable] = 
  {v: Variable _ NEW[VariableRec];
   m _ m + 1;
   v.i _ 0;
   v.j _ m;
   x[m] _ v;
   FOR i: INT IN [1 .. n] DO matrix[i][m] _ 0 ENDLOOP;
   RETURN [v]};

Pivot: PUBLIC PROC[k, l: Variable] = 
  {k: INT _ rowvar.i;
    l: INT _ rowvar.j;
    -- We will pivot at the matrix element matrix[k, l].
    i, j: INT; -- These are temporaries in the loops below.
   {temp:Variable _ y[k]; y[k] _ x[l]; x[l] _ temp};  -- swap the row and column headers
   {pp: REAL _ 1.0 / matrix[k][l]; matrix[k][l] _ 1.0;
     FOR j IN [1..m] DO 
       matrix[k][j] _ matrix[k][j] * pp -- divide everything in pivot row by the pivot element:
     ENDLOOP;
     FOR i IN [1..n] DO 
       IF i # k THEN 
        FOR j IN [1..m] DO 
          IF j # l THEN 
          -- for each matrix[i,j] outside the pivot row and column:
          matrix[i][j] _ matrix[i][j] - matrix[i][l] * matrix[k][j]; 
            -- note that matrix[k,j] was already * pp.
       ENDLOOP 
     ENDLOOP;
     -- Finally process pivot column:
     FOR i IN [1..n] 
       DO IF i # k THEN matrix[i][l] _ -matrix[i][l] * pp 
     ENDLOOP}};

AssertZ: PROC[l: LinearForm, PropEQ: PROC[u, v: Variable]] =
  {IF NOT satisfiable THEN ERROR;  
   {temp: Variable _ NewRowVar[];
    WHILE l # NIL DO Add[to: temp, what: l.first]; l _ l.rest ENDLOOP;
    -- now temp is a row variable that represents the linear form l
    Kill[temp, PropEQ]}};

NewRowVar: PROC[] RETURNS [v: Variable] =
  {v _ NEW[VariableRec];
   n _ n + 1;
   v.i _ n;
   v.j _ 0;
   y[n] _ v;
   FOR j: INT IN [1 .. m] DO matrix[n][j] _ 0 ENDLOOP;
   RETURN[v]};

Add: PROC[to: Variable, what: LinearMonomial] =
 -- to is assumed to be a row variable.
 {IF what.var.i = 0 
  THEN matrix[to.i][what.var.j] _ matrix[to.i][what.var.j] + what.coef
  ELSE FOR j: INT IN [1 .. m] 
           DO matrix[to.i][j] _ matrix[to.i][j] + what.coef * matrix[what.var.i][j]
           ENDLOOP};

Kill: PROC[v: Variable, PropEQ: PROC[u, v: Variable]] =
  {u, champion: Variable;
   pec: PivotEffectCode;
   IF matrix[v.i][unit.j] < 0 
   THEN FOR j: INT IN [1 .. m] DO matrix[v.i][j] _ - matrix[v.i][j] ENDLOOP;
     -- Asserting -v = 0 is the same as asserting v = 0.
     -- Now the sample value of v is non-positive.  Thus maximizing v will
     -- move it towards zero.
   [u, champion, pec] _ MaximalPivot[v];
   SELECT TRUE FROM
      pec = ToPivotAtUWouldMakeVPositive OR
      pec = VIsManifestlyUnbounded =>  
         {Pivot[champion.i, u.j]; KillColVar[champion]; CheckRowEqualities[PropEQ]};
      pec = VIsMaximizedAtZero => {KillAllCols[v, PropEQ]; CheckRowEqualities[PropEQ]};
      pec = VIsMaximizedAndNegative => satisfiable _ FALSE;
   ENDCASE => ERROR};

KillColVar: PROC[v: Variable, PropEQ: PROC[u, v: Variable]] =
  {IF v.j # m 
   THEN {FOR i: INT IN [1 .. n] DO matrix[i][v.j] _ matrix[i][m] ENDLOOP;
            x[v.j] _ x[m]};
   m _ m - 1};

KillAllCols: PROC[v:Variable, PropEQ: PROC[u, v: Variable]] =
 -- kill all col vars that have a non-zero entry in v's row.
 {j: INT _ 1;
  WHILE j <= m DO
    IF j # unit.j THEN KillColVar[x[j], PropEQ];
    j _ j + 1
  ENDLOOP};

PivotEffectCode: TYPE = 
  {ToPivotAtUWouldMakeVPositive, VIsMaximizedAtZero, VIsMaximizedAndNegative,
   ToPivotAtUWouldIncreaseVButNotMakeItPositive, VIsManifestlyUnbounded};
   -- The last code is misleadingly named: it may "increase" v by zero.
   -- To guarantee termination, some cycle-avoding 
 
MaximalPivot: PROC[v: Variable] RETURNS [u: Variable, champion: Variable _ NIL, pec: PivotEffectCode] =
  {[u, champion, pec] _ FindPivot[v];
   WHILE pec = ToPivotAtUWouldIncreaseVButNotMakeItPositive
   DO Pivot[champion.i, u.j]; [u, champion, pec] _ FindPivot[v] ENDLOOP};

FindPivot: PROC[v: Variable] RETURNS [u: Variable, champion: Variable _ NIL, pec: PivotEffectCode] =
  {j: INT _ 0;
   sgn: REAL;
   FOR jj: INT IN [1 .. m] 
   DO IF (matrix[v.i][jj] # 0 AND NOT y[i].restricted) OR matrix[v.i][jj] > 0 
       THEN {j _ jj; sgn _ - SGN[matrix[v.i][jj]]}
   ENDLOOP;   
   IF j = 0 THEN 
     SELECT TRUE FROM
       matrix[v.i][unit.j] < 0 => {pec _ VIsMaximizedAndNegative; RETURN};
       matrix[v.i][unit.j] = 0 => {pec _ VIsMaximizedAtZero; RETURN}
     ENDCASE => ERROR;
   u _ x[j]; -- we will pivot in column j
  {champ: INT;
   score: REAL;
   champ _ v.i;
   score _ Real.PlusInfinity;
   FOR ii: INT IN [1..n] 
   DO IF ii # v.i AND y[ii].restricted AND sgn * matrix[ii][j] < 0
       THEN IF ABS[matrix[ii][unit.j] / matrix[ii][j]] < score
               THEN {champ _ ii; score _ ABS[matrix[ii][unit.j] / matrix[ii][j]]}
   ENDLOOP;
   champion _ y[champ];
   IF champ = i THEN {pec _ VIsManifestlyUnbounded; RETURN};
   {newSampleValueOfV: REAL = 
     matrix[v.i][unit.j] - matrix[v.i][j] * matrix[i][unit.j] / matrix[i][j];
    SELECT TRUE FROM
      newSampleValueOfV <= 0 => pec _ ToPivotAtUWouldIncreaseVButNotMakeItPositive;
      newSampleValueOfV > 0  => pec _ ToPivotAtUWouldMakeVPositive
    ENDCASE => ERROR}}};

Restrict: PROC[v:Variable, PropEQ: PROC[u, v: Variable]] =
  {IF v.restricted THEN RETURN;
   MakeVOwnARow[v];
   -- Now v owns a row, or else it owns an empty column.
   IF v.i = 0 THEN {v.restricted _ TRUE; RETURN};
   IF matrix[v.i][unit.j] < 0 THEN {v.restricted _ TRUE; RETURN};
   {u: Variable; pec: PivotEffectCode;
    [u, champion, pec] _ MaximalPivot[v];
    SELECT TRUE FROM
      pec = ToPivotAtUWouldMakeVPositive => {v.restricted _ TRUE; Pivot[champion.i, u.j]};
      pec = VIsMaximizedAtZero => {KillAllCols[champion, PropEQ]};
      pec = VIsMaximizedAndNegative => {satisfiable _ FALSE};
      pec = VIsManifestlyUnbounded => {v.restricted _ TRUE; Pivot[v.i, u.j]};
    ENDCASE => ERROR}};
  
 Maximize: PROC[v:Variable] RETURNS [unbounded: BOOLEAN _ FALSE] = {
 u, champion: Variable;
 pec: PivotEffectCode;
 IF v.i # 0 THEN MakeVOwnARow[v]; -- try to make v a row
 IF v.i = 0 THEN RETURN[NOT v.restricted]; 
 [u, champion, pec] _ FindPivot[v];
 DO SELECT TRUE FROM
      pec = ToPivotAtUWouldMakeVPositive OR
      pec = ToPivotAtUWouldIncreaseVButNotMakeItPositive => {
          Pivot[champion.i, u.j]; [u, champion, pec] _ FindPivot[v]};
    ENDCASE => EXIT; ENDLOOP;
 unbounded _ pec = VIsManifestlyUnbounded;
 };
 
 MakeVOwnARow: PROC [v:Variable] = {
   champ: INT;
   score: REAL;
   sgn: REAL;
   i: INT _ 1;
   WHILE i < n AND matrix[i][v.j] = 0 DO i _ i + 1 ENDLOOP;
   SELECT TRUE FROM
     i = n => RETURN;
     matrix[i][v.j] > 0 => sgn _ +1;
     matrix[i][v.j] < 0 => sgn _ -1
   ENDCASE => ERROR;
   champ _ v.i;
   score _ Real.PlusInfinity;
   FOR ii: INT IN [1..n] 
   DO IF ii # v.i AND y[ii].restricted AND sgn * matrix[ii][j] < 0
       THEN IF ABS[matrix[ii][unit.j] / matrix[ii][j]] < score
               THEN {champ _ ii; score _ ABS[matrix[ii][unit.j] / matrix[ii][j]]}
   ENDLOOP;
   champion _ y[champ];
   Pivot[champ.i, v.j];
   };
  
 }.