/*****
* simplex.asy
* Andy Hammerlindl 2004/07/27
*
* Solves the two-variable linear programming problem using the simplex method.
* This problem is specialized in that the second variable, "b", does not have
* a non-negativity condition, and the first variable, "a", is the quantity
* being maximized.
* Correct execution of the algorithm also assumes that the coefficient of "b"
* will be +1 or -1 in every added restriction, and that the problem can be
* initialized to a valid state by pivoting b with one of the slack
* variables.  This assumption may in fact be incorrect.
*****/

struct problem {
 typedef int var;
 static var VAR_A = 0;
 static var VAR_B = 1;

 static int OPTIMAL = -1;
 static var UNBOUNDED = -2;
 static int INVALID = -3;

 struct row {
   real c, t[];
 }

 // The variables of the rows.
 // Initialized for the two variable problem.
 var[] v = {VAR_A, VAR_B};

 // The rows of equalities.
 row rowA() {
   row r = new row;
   r.c = 0;
   r.t = new real[] {1, 0};
   return r;
 }
 row rowB() {
   row r = new row;
   r.c = 0;
   r.t = new real[] {0, 1};
   return r;
 }
 row[] rows = {rowA(), rowB()};

 // The number of original variables.
 int n = rows.length;

 // Pivot the variable v[col] with vp.
 void pivot(int col, var vp)
 {
   var vc=v[col];

   // Recalculate rows v[col] and vp for the pivot-swap.
   row rvc = rows[vc], rvp = rows[vp];
   real factor=1/rvp.t[col]; // NOTE: Handle rvp.t[col] == 0 case.
   rvc.c=-rvp.c*factor;
   rvp.c=0;
   rvc.t=-rvp.t*factor;
   rvp.t *= 0;
   rvc.t[col]=factor;
   rvp.t[col]=1;

   var a=min(vc,vp);
   var b=max(vc,vp);

   // Recalculate the rows other than the two used for the above pivot.
   for (var i = 0; i < a; ++i) {
     row r=rows[i];
     real m = r.t[col];
     r.c += m*rvc.c;
     r.t += m*rvc.t;
     r.t[col]=m*factor;
   }
   for (var i = a+1; i < b; ++i) {
     row r=rows[i];
     real m = r.t[col];
     r.c += m*rvc.c;
     r.t += m*rvc.t;
     r.t[col]=m*factor;
   }
   for (var i = b+1; i < rows.length; ++i) {
     row r=rows[i];
     real m = r.t[col];
     r.c += m*rvc.c;
     r.t += m*rvc.t;
     r.t[col]=m*factor;
   }

   // Relabel the vars.
   v[col] = vp;
 }

 // As b does not have a non-negativity condition, it must initially be
 // pivoted out for a variable that does.  This selects the initial
 // variable to pivot with b.  It also assumes that there is a valid
 // solution with a == 0 to the linear programming problem, and if so, it
 // picks a pivot to get to that state.  In our case, a == 0 corresponds to
 // a picture with the user coordinates shrunk down to zero, and if that
 // doesn't fit, nothing will.
 //
 // If b has a minimal value, choose a pivot that will give b its minimal
 // value.  Otherwise, if b has maximal value, choose a pivot to give b its
 // maximal value.
 var initVar()
 {
   real min;
   var argmin;
   var i=2;
   for (; i < rows.length; ++i) {
     row r=rows[i];
     if (r.t[VAR_B] > 0) {
       min=r.c/r.t[VAR_B];
       argmin=i;
       break;
     }
   }
   for (; i < rows.length; ++i) {
     row r=rows[i];
     if (r.t[VAR_B] > 0) {
       real val=r.c/r.t[VAR_B];
       if (val < min) {
         min=val;
         argmin=i;
       }
     }
   }

   if(argmin != 0) return argmin;

   real max;
   var argmax;
   var i=2;
   for (; i < rows.length; ++i) {
     row r=rows[i];
     if (r.t[VAR_B] < 0) {
       max=r.c/r.t[VAR_B];
       argmax=i;
       break;
     }
   }

   for (; i < rows.length; ++i) {
     row r=rows[i];
     if (r.t[VAR_B] < 0) {
       real val=r.c/r.t[VAR_B];
       if (val > max) {
         max=val;
         argmax=i;
       }
     }
   }

   if(argmax != 0) return argmax;

   return UNBOUNDED;
 }

 // Initialize the linear program problem by moving into an acceptable state
 // this assumes that b is unrestrained and is the second variable.
 // NOTE: Works in limited cases, may be bug-ridden.
 void init()
 {
   // Find the lowest constant term in the equations.
   var lowest = 0;
   for (var i = 2; i < rows.length; ++i) {
     if (rows[i].c < rows[lowest].c)
       lowest = i;
   }

   // Pivot if necessary.
   if (lowest != 0)
     pivot(VAR_B, lowest);
 }

 // Selects a column to pivot on.  Returns OPTIMAL if the current state is
 // optimal.  Assumes we are optimizing the first row.
 int selectColumn()
 {
   int i=find(rows[0].t > 0,1);
   return (i >= 0) ? i : OPTIMAL;
 }

 // Select the new variable associated with a pivot on the column given.
 // Returns UNBOUNDED if the space is unbounded.
 var selectVar(int col)
 {
   // We assume that the first two vars (a and b) once swapped out, won't be
   // swapped back in.  This finds the variable which gives the tightest
   // non-negativity condition restricting our optimization.  This turns
   // out to be the max of c/t[col].  Note that as c is positive, and
   // t[col] is negative, all c/t[col] will be negative, so we are finding
   // the smallest in magnitude.
   var vp=UNBOUNDED;
   real max=0;
   int i=2;
   for (; i < rows.length; ++i) {
     row r=rows[i];
     if(r.c < 0) r.c=0; // Fix any numerical precision error
     if(r.t[col] < 0) {
       max=r.c/r.t[col]; vp=i;
       break;
     }
   }
   for (; i < rows.length; ++i) {
     row r=rows[i];
     if(r.c < 0) r.c=0; // Fix any numerical precision error
     if(r.c < max*r.t[col]) {
       max=r.c/r.t[col]; vp=i;
     }
   }

   return vp;
 }

 // Checks that the rows are in a valid state.
 bool valid()
 {
   // Checks that constants are valid.
   bool validConstants() {
     for (int i = 0; i < rows.length; ++i)
       // Do not test the row for b, as it does not have a non-negativity
       // condition.
       if (i != VAR_B && rows[i].c < 0)
          return false;
     return true;
   }

   // Check a variable to see if its row is simple.
   // NOTE: Simple rows could be optimized out, since they are not really
   // used.
   bool validVar(int col) {

     var vc = v[col];
     row rvc = rows[vc];

     if (rvc.c != 0)
       return false;
     for (int i = 0; i < n; ++i)
       if (rvc.t[i] != (i == col ? 1 : 0))
         return false;

     return true;
   }

   if (!validConstants()) {
     return false;
   }
   for (int i = 0; i < n; ++i)
     if (!validVar(i)) {
       return false;
     }

   return true;
 }


 // Perform the algorithm to find the optimal solution.  Returns OPTIMAL,
 // UNBOUNDED, or INVALID (if no solution is possible).
 int optimize()
 {
   // Put into a valid state to begin and pivot b out.
   var iv=initVar();
   if (iv == UNBOUNDED)
     return iv;
   pivot(VAR_B, iv);

   if (!valid())
     return INVALID;

   while(true) {
     int col = selectColumn();

     if (col == OPTIMAL)
       return col;
     var vp = selectVar(col);

     if (vp == UNBOUNDED)
       return vp;

     pivot(col, vp);
   }

   // Shouldn't reach here.
   return INVALID;
 }

 // Add a restriction to the problem:
 // t1*a + t2*b + c >= 0
 void addRestriction(real t1, real t2, real c)
 {
   row r = new row;
   r.c = c;
   r.t = new real[] {t1, t2};
   rows.push(r);
 }

 // Return the value of a computed.
 real a()
 {
   return rows[VAR_A].c;
 }

 // Return the value of b computed.
 real b()
 {
   return rows[VAR_B].c;
 }
}

You are viewing proxied material from ftp.icm.edu.pl. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.