// Rational simplex solver written by John C. Bowman and Pouria Ramazi, 2018.
import rational;

bool optimizeTableau=true;

int[] artificialColumn;

void simplexInit(rational[] c, rational[][] A, int[] s=new int[],
                rational[] b, int count) {}
void simplexTableau(rational[][] E, int[] Bindices, int I=-1, int J=-1,
                   int n=E[0].length-1) {}
void simplexPhase1(rational[] c, rational[][] A, rational[] b,
                  int[] Bindices) {}
void simplexPhase2() {}

void simplexWrite(rational[][] E, int[] Bindices, int, int)
{
 int m=E.length-1;
 int n=E[0].length-1;

 write(E[m][0],tab);
 for(int j=1; j <= n; ++j)
   write(E[m][j],tab);
 write();

 for(int i=0; i < m; ++i) {
   write(E[i][0],tab);
   for(int j=1; j <= n; ++j) {
     write(E[i][j],tab);
   }
   write();
 }
 write();
};

struct simplex {
 static int OPTIMAL=0;
 static int UNBOUNDED=1;
 static int INFEASIBLE=2;

 int case;
 rational[] x;
 rational[] xStandard;
 rational cost;
 rational[] d;
 bool dual=false;

 int m,n;
 int J;

 // Row reduce based on pivot E[I][J]
 void rowreduce(rational[][] E, int N, int I, int J) {
   rational[] EI=E[I];
   rational v=EI[J];
   for(int j=0; j < J; ++j) EI[j] /= v;
   EI[J]=1;
   for(int j=J+1; j <= N; ++j) EI[j] /= v;

   for(int i=0; i < I; ++i) {
     rational[] Ei=E[i];
     rational EiJ=Ei[J];
     for(int j=0; j < J; ++j)
       Ei[j] -= EI[j]*EiJ;
     Ei[J]=0;
     for(int j=J+1; j <= N; ++j)
       Ei[j] -= EI[j]*EiJ;
   }
   for(int i=I+1; i <= m; ++i) {
     rational[] Ei=E[i];
     rational EiJ=Ei[J];
     for(int j=0; j < J; ++j)
       Ei[j] -= EI[j]*EiJ;
     Ei[J]=0;
     for(int j=J+1; j <= N; ++j)
       Ei[j] -= EI[j]*EiJ;
   }
 }

 int iterate(rational[][] E, int N, int[] Bindices, bool phase1=false) {
   while(true) {
     // Bland's rule: first negative entry in reduced cost (bottom) row enters
     rational[] Em=E[m];
     for(J=1; J <= N; ++J)
       if(Em[J] < 0) break;

     if(J > N)
       break;

     int I=-1;
     rational t;
     for(int i=0; i < m; ++i) {
       rational u=E[i][J];
       if(u > 0) {
         t=E[i][0]/u;
         I=i;
         break;
       }
     }
     for(int i=I+1; i < m; ++i) {
       rational u=E[i][J];
       if(u > 0) {
         rational r=E[i][0]/u;
         if(r <= t && (r < t || Bindices[i] < Bindices[I])) {
           t=r; I=i;
         } // Bland's rule: exiting variable has smallest minimizing subscript
       }
     }
     if(I == -1)
       return UNBOUNDED; // Can only happen in Phase 2.

     simplexTableau(E,Bindices,I,J);

     // Generate new tableau
     Bindices[I]=J;
     rowreduce(E,N,I,J);

     if(phase1 && Em[0] == 0) break;
   }
   return OPTIMAL;
 }

 int iterateDual(rational[][] E, int N, int[] Bindices, bool phase1=false) {
   while(true) {
     // Bland's rule: negative variable with smallest subscript exits
     int I;
     for(I=0; I < m; ++I) {
       if(E[I][0] < 0) break;
     }

     if(I == m)
       break;

     for(int i=I+1; i < m; ++i) {
       if(E[i][0] < 0 && Bindices[i] < Bindices[I])
         I=i;
     }

     rational[] Em=E[m];
     rational[] EI=E[I];
     int J=0;
     rational t;
     for(int j=1; j <= N; ++j) {
       rational u=EI[j];
       if(u < 0) {
         t=-Em[j]/u;
         J=j;
         break;
       }
     }
     for(int j=J+1; j <= N; ++j) {
       rational u=EI[j];
       if(u < 0) {
         rational r=-Em[j]/u;
         if(r <= t && (r < t || j < J)) {
           t=r; J=j;
         } // Bland's rule: smallest minimizing subscript enters
       }
     }
     if(J == 0)
       return INFEASIBLE; // Can only happen in Phase 2.

     simplexTableau(E,Bindices,I,J);

     // Generate new tableau
     Bindices[I]=J;
     rowreduce(E,N,I,J);
   }
   return OPTIMAL;
 }

 // Try to find a solution x to Ax=b that minimizes the cost c^T x,
 // where A is an m x n matrix, x is a vector of n non-negative numbers,
 // b is a vector of length m, and c is a vector of length n.
 // Can set phase1=false if the last m columns of A form the identity matrix.
 void operator init(rational[] c, rational[][] A, rational[] b,
                    bool phase1=true) {
   // Phase 1
   m=A.length;
   if(m == 0) {case=INFEASIBLE; return;}
   n=A[0].length;
   if(n == 0) {case=INFEASIBLE; return;}

   rational[][] E=new rational[m+1][n+1];
   rational[] Em=E[m];

   for(int j=1; j <= n; ++j)
     Em[j]=0;

   for(int i=0; i < m; ++i) {
     rational[] Ai=A[i];
     rational[] Ei=E[i];
     if(b[i] >= 0 || dual) {
       for(int j=1; j <= n; ++j) {
         rational Aij=Ai[j-1];
         Ei[j]=Aij;
         Em[j] -= Aij;
       }
     } else {
       for(int j=1; j <= n; ++j) {
         rational Aij=-Ai[j-1];
         Ei[j]=Aij;
         Em[j] -= Aij;
       }
     }
   }

   void basicValues() {
     rational sum=0;
     for(int i=0; i < m; ++i) {
       rational B=dual ? b[i] : abs(b[i]);
       E[i][0]=B;
       sum -= B;
     }
     Em[0]=sum;
   }

   int[] Bindices;

   if(phase1) {
     Bindices=new int[m];
     int k=0;

     artificialColumn.delete();
     // Check for redundant basis vectors.
     for(int p=0; p < m; ++p) {
       bool checkBasis(int j) {
         for(int i=0; i < m; ++i) {
           rational[] Ei=E[i];
           if(i != p ? Ei[j] != 0 : Ei[j] <= 0)
             return false;
         }
         return true;
       }

       int checkTableau() {
         if(optimizeTableau)
           for(int j=1; j <= n; ++j)
             if(checkBasis(j)) return j;
         return 0;
       }

       int j=checkTableau();
       Bindices[p]=n+1+p;
       if(j == 0) { // Add an artificial variable
         artificialColumn.push(p+1);
         for(int i=0; i < p; ++i)
           E[i].push(0);
         E[p].push(1);
         for(int i=p+1; i < m; ++i)
           E[i].push(0);
         E[m].push(0);
         ++k;
       }
     }

     basicValues();

     simplexPhase1(c,A,b,Bindices);

     iterate(E,n+k,Bindices,true);

     if(Em[0] != 0) {
       simplexTableau(E,Bindices);
     case=INFEASIBLE;
     return;
     }
   } else {
      Bindices=sequence(new int(int x){return x;},m)+n-m+1;
      basicValues();
   }

   rational[] cB=phase1 ? new rational[m] : c[n-m:n];
   rational[][] D=phase1 ? new rational[m+1][n+1] : E;
   if(phase1) {
     // Drive artificial variables out of basis.
     for(int i=0; i < m; ++i) {
       if(Bindices[i] > n) {
         rational[] Ei=E[i];
         int j;
         for(j=1; j <= n; ++j)
           if(Ei[j] != 0) break;
         if(j > n) continue;
         simplexTableau(E,Bindices,i,j,n);
         Bindices[i]=j;
         rowreduce(E,n,i,j);
       }
     }
     simplexTableau(E,Bindices,-1,-1,n);
     int ip=0; // reduced i
     for(int i=0; i < m; ++i) {
       int k=Bindices[i];
       if(k > n) continue;
       Bindices[ip]=k;
       cB[ip]=c[k-1];
       rational[] Dip=D[ip];
       rational[] Ei=E[i];
       for(int j=1; j <= n; ++j)
         Dip[j]=Ei[j];
       Dip[0]=Ei[0];
       ++ip;
     }

     rational[] Dip=D[ip];
     rational[] Em=E[m];
     for(int j=1; j <= n; ++j)
       Dip[j]=Em[j];
     Dip[0]=Em[0];

     if(m > ip) {
       Bindices.delete(ip,m-1);
       D.delete(ip,m-1);
       m=ip;
     }
   }

   rational[] Dm=D[m];
   for(int j=1; j <= n; ++j) {
     rational sum=0;
     for(int k=0; k < m; ++k)
       sum += cB[k]*D[k][j];
     Dm[j]=c[j-1]-sum;
   }

   rational sum=0;
   for(int k=0; k < m; ++k)
     sum += cB[k]*D[k][0];
   Dm[0]=-sum;

   simplexPhase2();

   case=(dual ? iterateDual : iterate)(D,n,Bindices);
   simplexTableau(D,Bindices);

   if(case != INFEASIBLE) {
     x=new rational[n];
     for(int j=0; j < n; ++j)
       x[j]=0;

     for(int k=0; k < m; ++k)
       x[Bindices[k]-1]=D[k][0];

     xStandard=copy(x);
   }

   if(case == UNBOUNDED) {
     d=new rational[n];
     for(int j=0; j < n; ++j)
       d[j]=0;
     d[J-1]=1;
     for(int k=0; k < m; ++k)
       d[Bindices[k]-1]=-D[k][J];
   }

   if(case != OPTIMAL)
     return;

   cost=-Dm[0];
 }

 // Try to find a solution x to sgn(Ax-b)=sgn(s) that minimizes the cost
 // c^T x, where A is an m x n matrix, x is a vector of n non-negative
 // numbers, b is a vector of length m, and c is a vector of length n.
 void operator init(rational[] c, rational[][] A, int[] s, rational[] b) {
   int m=A.length;
   if(m == 0) {case=INFEASIBLE; return;}
   int n=A[0].length;
   if(n == 0) {case=INFEASIBLE; return;}

   int count=0;
   for(int i=0; i < m; ++i)
     if(s[i] != 0) ++count;

   rational[][] a=new rational[m][n+count];

   for(int i=0; i < m; ++i) {
     rational[] ai=a[i];
     rational[] Ai=A[i];
     for(int j=0; j < n; ++j) {
       ai[j]=Ai[j];
     }
   }

   int k=0;

   bool phase1=false;
   dual=count == m && all(c >= 0);

   for(int i=0; i < m; ++i) {
     rational[] ai=a[i];
     for(int j=0; j < k; ++j)
       ai[n+j]=0;
     int si=s[i];
     if(k < count)
       ai[n+k]=-si;
     for(int j=k+1; j < count; ++j)
       ai[n+j]=0;
     if(si == 0) phase1=true;
     else {
       ++k;
       rational bi=b[i];
       if(bi == 0) {
         if(si == 1) {
           s[i]=-1;
           for(int j=0; j < n+count; ++j)
             ai[j]=-ai[j];
         }
       } else if(dual && si == 1) {
         b[i]=-bi;
         s[i]=-1;
         for(int j=0; j < n+count; ++j)
           ai[j]=-ai[j];
       } else if(si*bi > 0)
         phase1=true;
     }
   }

   if(dual) phase1=false;
   rational[] C=concat(c,array(count,rational(0)));
   simplexInit(C,a,b,count);
   operator init(C,a,b,phase1);

   if(case != INFEASIBLE && count > 0)
     x.delete(n,n+count-1);
 }
}

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.