option randseed 0;

# primal simplex method


param m;
param n;
param bb {1..m};
param cc {1..n};
param AA {1..m, 1..n};
param A {1..m, 1..n};
param b {1..m};
param c {1..n};
param nonbasics {1..n};
param basics {1..m};

param xx;
param yy;
param obj;

param iter;
param maxc;
param maxaoverb;
param col;
param row;
param Arow {1..n};
param Acol {1..m};
param a;
param brow;
param ccol;
param ii;
param jj;

let m := 10;
let n := 2;

param x0 := 4*Uniform01();
param y0 := 4*Uniform01();
param r := 5;
param pi := 4*atan(1);
param theta {j in 1..m} := (3.*Uniform01()-1.)*pi/2.;
param x1 {j in 1..m} := x0 + r*cos(theta[j]);
param y1 {j in 1..m} := y0 + r*sin(theta[j]);
let {i in 1..m} A[i,1] := x1[i] - x0;
let {i in 1..m} A[i,2] := y1[i] - y0;
let {i in 1..m} b[i] := (x1[i]-x0)*x1[i] + (y1[i]-y0)*y1[i];
let {i in 1..m, j in 1..n} A[i,j] := ceil(20*A[i,j]);
let {i in 1..m} b[i] := round(20*b[i]);
let {j in 1..n} c[j] := ceil(9*Uniform01());

# save the original A, b, c data
let {i in 1..m, j in 1..n} AA[i,j] := A[i,j];
let {i in 1..m} bb[i] := b[i];
let {j in 1..n} cc[j] := c[j];

let {i in 1..m, j in 1..n} A[i,j] := -A[i,j];
let {j in 1..n} nonbasics[j] := j;
let {i in 1..m} basics[i] := n+i;

let iter := 0;
let obj := 0;
printf "(x,y) = %8f, %8f,   objective function = %8f \n", 0,0,0;
repeat while (max {j in 1..n} c[j]) > 0 {
   let maxc := 0;
   for {j in 1..n} {
	  if (c[j] > maxc) then {
	     let maxc := c[j];
	     let col := j;
	  }
   }
   let maxaoverb := 0;
   for {i in 1..m} {
	  if (-A[i,col]/b[i] > maxaoverb) then {
	     let maxaoverb := -A[i,col]/b[i];
	     let row := i;
	  }
   }

   let {j in 1..n} Arow[j] := A[row,j];
   let {i in 1..m} Acol[i] := A[i,col];
   let a := A[row,col];
   let {i in 1..m, j in 1..n} A[i,j] := A[i,j] - Acol[i]*Arow[j]/a;
   let {j in 1..n} A[row,j] := -Arow[j]/a;
   let {i in 1..m} A[i,col] :=  Acol[i]/a;
   let A[row,col] := 1/a;

   let brow := b[row];
   let {i in 1..m} b[i] := b[i] - brow*Acol[i]/a;
   let b[row] := -brow/a;

   let ccol := c[col];
   let {j in 1..n} c[j] := c[j] - ccol*Arow[j]/a;
   let c[col] := ccol/a;

   let jj := nonbasics[col];
   let ii := basics[row];
   let basics[row] := jj;
   let nonbasics[col] := ii;

   let iter := iter+1;

   let obj := obj - ccol*brow/a;
   let xx := 0;
   let yy := 0;
   for {i in 1..m} {
           if (basics[i] == 1) then {let xx := b[i];}
           if (basics[i] == 2) then {let yy := b[i];}
   }
   printf "(x,y) = %8f, %8f,   objective function = %8f \n", xx,yy, obj;
}
printf "  (m,n) = %3d, %3d,   iters = %3d \n", m, n, iter;