// Real simplex solver written by John C. Bowman and Pouria Ramazi, 2018. struct simplex { static int OPTIMAL=0; static int UNBOUNDED=1; static int INFEASIBLE=2; int case; real[] x; real cost; bool dual=false; int m,n; int J; real EpsilonA; // Row reduce based on pivot E[I][J] void rowreduce(real[][] E, int N, int I, int J) { real[] EI=E[I]; real v=EI[J]; for(int j=0; j < J; ++j) EI[j] /= v; EI[J]=1.0; for(int j=J+1; j <= N; ++j) EI[j] /= v; for(int i=0; i < I; ++i) { real[] Ei=E[i]; real EiJ=Ei[J]; for(int j=0; j < J; ++j) Ei[j] -= EI[j]*EiJ; Ei[J]=0.0; for(int j=J+1; j <= N; ++j) Ei[j] -= EI[j]*EiJ; } for(int i=I+1; i <= m; ++i) { real[] Ei=E[i]; real EiJ=Ei[J]; for(int j=0; j < J; ++j) Ei[j] -= EI[j]*EiJ; Ei[J]=0.0; for(int j=J+1; j <= N; ++j) Ei[j] -= EI[j]*EiJ; } } int iterate(real[][] E, int N, int[] Bindices, bool phase1=false) { while(true) { // Bland's rule: first negative entry in reduced cost (bottom) row enters real[] Em=E[m]; for(J=1; J <= N; ++J) if(Em[J] < 0) break; if(J > N) break; int I=-1; real t; for(int i=0; i < m; ++i) { real u=E[i][J]; if(u > EpsilonA) { t=E[i][0]/u; I=i; break; } } for(int i=I+1; i < m; ++i) { real u=E[i][J]; if(u > EpsilonA) { real 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. // Generate new tableau Bindices[I]=J; rowreduce(E,N,I,J); if(phase1 && abs(Em[0]) <= EpsilonA) break; } return OPTIMAL; } int iterateDual(real[][] 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; } real[] Em=E[m]; real[] EI=E[I]; int J=0; real t; for(int j=1; j <= N; ++j) { real u=EI[j]; if(u < -EpsilonA) { t=-Em[j]/u; J=j; break; } } for(int j=J+1; j <= N; ++j) { real u=EI[j]; if(u < -EpsilonA) { real r=-Em[j]/u; if(r < t) { t=r; J=j; } // Bland's rule: smallest minimizing subscript enters } } if(J == 0) return INFEASIBLE; // Can only happen in Phase 2. // 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(real[] c, real[][] A, real[] b, bool phase1=true) { static real epsilon=sqrt(realEpsilon); real normA=norm(A); real epsilonA=100.0*realEpsilon*normA; EpsilonA=epsilon*normA; // Phase 1 m=A.length; if(m == 0) {case=INFEASIBLE; return;} n=A[0].length; if(n == 0) {case=INFEASIBLE; return;} real[][] E=new real[m+1][n+1]; real[] Em=E[m]; for(int j=1; j <= n; ++j) Em[j]=0; for(int i=0; i < m; ++i) { real[] Ai=A[i]; real[] Ei=E[i]; if(b[i] >= 0 || dual) { for(int j=1; j <= n; ++j) { real Aij=Ai[j-1]; Ei[j]=Aij; Em[j] -= Aij; } } else { for(int j=1; j <= n; ++j) { real Aij=-Ai[j-1]; Ei[j]=Aij; Em[j] -= Aij; } } } void basicValues() { real sum=0; for(int i=0; i < m; ++i) { real 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; // Check for redundant basis vectors. for(int p=0; p < m; ++p) { bool checkBasis(int j) { for(int i=0; i < m; ++i) { real[] Ei=E[i]; if(i != p ? abs(Ei[j]) >= epsilonA : Ei[j] <= epsilonA) return false; } return true; } int checkTableau() { 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 for(int i=0; i < p; ++i) E[i].push(0.0); E[p].push(1.0); for(int i=p+1; i < m; ++i) E[i].push(0.0); E[m].push(0.0); ++k; } } basicValues(); iterate(E,n+k,Bindices,true); if(abs(Em[0]) > EpsilonA) { case=INFEASIBLE; return; } } else { Bindices=sequence(new int(int x){return x;},m)+n-m+1; basicValues(); } real[] cB=phase1 ? new real[m] : c[n-m:n]; real[][] D=phase1 ? new real[m+1][n+1] : E; if(phase1) { // Drive artificial variables out of basis. for(int i=0; i < m; ++i) { if(Bindices[i] > n) { real[] Ei=E[i]; int j; for(j=1; j <= n; ++j) if(abs(Ei[j]) > EpsilonA) break; if(j > n) continue; Bindices[i]=j; rowreduce(E,n,i,j); } } 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]; real[] Dip=D[ip]; real[] Ei=E[i]; for(int j=1; j <= n; ++j) Dip[j]=Ei[j]; Dip[0]=Ei[0]; ++ip; } real[] Dip=D[ip]; real[] 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; } } real[] Dm=D[m]; for(int j=1; j <= n; ++j) { real sum=0; for(int k=0; k < m; ++k) sum += cB[k]*D[k][j]; Dm[j]=c[j-1]-sum; } real sum=0; for(int k=0; k < m; ++k) sum += cB[k]*D[k][0]; Dm[0]=-sum; case=(dual ? iterateDual : iterate)(D,n,Bindices); if(case != OPTIMAL) return; for(int j=0; j < n; ++j) x[j]=0; for(int k=0; k < m; ++k) x[Bindices[k]-1]=D[k][0]; 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(real[] c, real[][] A, int[] s, real[] 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; real[][] a=new real[m][n+count]; for(int i=0; i < m; ++i) { real[] ai=a[i]; real[] Ai=A[i]; for(int j=0; j < n; ++j) { ai[j]=Ai[j]; } } int k=0; bool phase1=false; bool dual=count == m && all(c >= 0); for(int i=0; i < m; ++i) { real[] 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; real 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; operator init(concat(c,array(count,0.0)),a,b,phase1); if(case == OPTIMAL && count > 0) x.delete(n,n+count-1); } }