// 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) {
      write("n=",n);
      write(Bindices);
      // 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);
  }
}