diff options
Diffstat (limited to 'Build/source/utils/asymptote/base/simplex.asy')
-rw-r--r-- | Build/source/utils/asymptote/base/simplex.asy | 363 |
1 files changed, 363 insertions, 0 deletions
diff --git a/Build/source/utils/asymptote/base/simplex.asy b/Build/source/utils/asymptote/base/simplex.asy new file mode 100644 index 00000000000..e6d0410b23c --- /dev/null +++ b/Build/source/utils/asymptote/base/simplex.asy @@ -0,0 +1,363 @@ +// 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) { + 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); + } + return OPTIMAL; + } + + int iterateDual(real[][] E, int N, int[] Bindices) { + 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 p=0; + + // Check for redundant basis vectors. + 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 k=0; + while(p < m) { + int j=checkTableau(); + if(j > 0) + Bindices[p]=j; + else { // Add an artificial variable + Bindices[p]=n+1+k; + 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; + } + ++p; + } + + basicValues(); + iterate(E,n+k,Bindices); + + 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) { + int k=Bindices[i]; + if(k > 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); + } +} |