diff options
Diffstat (limited to 'Build/source/utils/asymptote/base/simplex.asy')
-rw-r--r-- | Build/source/utils/asymptote/base/simplex.asy | 156 |
1 files changed, 91 insertions, 65 deletions
diff --git a/Build/source/utils/asymptote/base/simplex.asy b/Build/source/utils/asymptote/base/simplex.asy index 040621dd904..b9ff172077b 100644 --- a/Build/source/utils/asymptote/base/simplex.asy +++ b/Build/source/utils/asymptote/base/simplex.asy @@ -11,7 +11,7 @@ struct simplex { int m,n; int J; - real epsilonA; + real EpsilonA; // Row reduce based on pivot E[I][J] void rowreduce(real[][] E, int N, int I, int J) { @@ -45,26 +45,26 @@ struct simplex { while(true) { // Find first negative entry in bottom (reduced cost) row real[] Em=E[m]; - for(J=0; J < N; ++J) + for(J=1; J <= N; ++J) if(Em[J] < 0) break; - if(J == N) + if(J > N) break; int I=-1; real M; for(int i=0; i < m; ++i) { real e=E[i][J]; - if(e > epsilonA) { - M=E[i][N]/e; + if(e > EpsilonA) { + M=E[i][0]/e; I=i; break; } } for(int i=I+1; i < m; ++i) { real e=E[i][J]; - if(e > epsilonA) { - real v=E[i][N]/e; + if(e > EpsilonA) { + real v=E[i][0]/e; if(v < M) {M=v; I=i;} // Bland's rule: choose smallest argmin } } @@ -80,34 +80,34 @@ struct simplex { int iterateDual(real[][] E, int N, int[] Bindices) { while(true) { - // Find first negative entry in right (basic variable) column + // Find first negative entry in zeroth (basic variable) column real[] Em=E[m]; int I; for(I=0; I < m; ++I) { - if(E[I][N] < 0) break; + if(E[I][0] < 0) break; } if(I == m) break; - int J=-1; + int J=0; real M; - for(int j=0; j < N; ++j) { + for(int j=1; j <= N; ++j) { real e=E[I][j]; - if(e < epsilonA) { + if(e < -EpsilonA) { M=-E[m][j]/e; J=j; break; } } - for(int j=J+1; j < N; ++j) { + for(int j=J+1; j <= N; ++j) { real e=E[I][j]; - if(e < epsilonA) { + if(e < -EpsilonA) { real v=-E[m][j]/e; if(v < M) {M=v; J=j;} // Bland's rule: choose smallest argmin } } - if(J == -1) + if(J == 0) return INFEASIBLE; // Can only happen in Phase 2. // Generate new tableau @@ -125,7 +125,9 @@ struct simplex { bool dual=false) { if(dual) phase1=false; static real epsilon=sqrt(realEpsilon); - epsilonA=epsilon*norm(A); + real normA=norm(A); + real epsilonA=100.0*realEpsilon*normA; + EpsilonA=epsilon*normA; // Phase 1 m=A.length; @@ -133,78 +135,103 @@ struct simplex { n=A[0].length; if(n == 0) {case=INFEASIBLE; return;} - int N=phase1 ? n+m : n; - real[][] E=new real[m+1][N+1]; + real[][] E=new real[m+1][n+1]; real[] Em=E[m]; - for(int j=0; j < n; ++j) + 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=0; j < n; ++j) { - real Aij=Ai[j]; + for(int j=1; j <= n; ++j) { + real Aij=Ai[j-1]; Ei[j]=Aij; Em[j] -= Aij; } } else { - for(int j=0; j < n; ++j) { - real Aij=-Ai[j]; + for(int j=1; j <= n; ++j) { + real Aij=-Ai[j-1]; Ei[j]=Aij; Em[j] -= Aij; } } } - if(phase1) { - for(int i=0; i < m; ++i) { - real[] Ei=E[i]; - for(int j=0; j < i; ++j) - Ei[n+j]=0.0; - Ei[n+i]=1.0; - for(int j=i+1; j < m; ++j) - Ei[n+j]=0.0; + 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; } - real sum=0; - for(int i=0; i < m; ++i) { - real B=dual ? b[i] : abs(b[i]); - E[i][N]=B; - sum -= B; - } - Em[N]=sum; - - if(phase1) - for(int j=0; j < m; ++j) - Em[n+j]=0.0; - int[] Bindices; if(phase1) { - Bindices=sequence(new int(int x){return x;},m)+n; - iterate(E,N,Bindices); + 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[J]) > epsilonA) { + if(abs(Em[0]) > EpsilonA) { case=INFEASIBLE; return; } - } else Bindices=sequence(new int(int x){return x;},m)+n-m; - + } 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) { + if(k > n) { real[] Ei=E[i]; int j; - for(j=0; j < n; ++j) - if(Ei[j] != 0) break; - if(j == n) continue; + for(j=1; j <= n; ++j) + if(abs(Ei[j]) > EpsilonA) break; + if(j > n) continue; Bindices[i]=j; rowreduce(E,n,i,j); } @@ -212,22 +239,22 @@ struct simplex { int ip=0; // reduced i for(int i=0; i < m; ++i) { int k=Bindices[i]; - if(k >= n) continue; + if(k > n) continue; Bindices[ip]=k; - cB[ip]=c[k]; + cB[ip]=c[k-1]; real[] Dip=D[ip]; real[] Ei=E[i]; - for(int j=0; j < n; ++j) + for(int j=1; j <= n; ++j) Dip[j]=Ei[j]; - Dip[n]=Ei[N]; + Dip[0]=Ei[0]; ++ip; } real[] Dip=D[ip]; real[] Em=E[m]; - for(int j=0; j < n; ++j) + for(int j=1; j <= n; ++j) Dip[j]=Em[j]; - Dip[n]=Em[N]; + Dip[0]=Em[0]; if(m > ip) { Bindices.delete(ip,m-1); @@ -237,17 +264,17 @@ struct simplex { } real[] Dm=D[m]; - for(int j=0; j < n; ++j) { + 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]-sum; + Dm[j]=c[j-1]-sum; } real sum=0; for(int k=0; k < m; ++k) - sum += cB[k]*D[k][n]; - Dm[n]=-sum; + sum += cB[k]*D[k][0]; + Dm[0]=-sum; case=(dual ? iterateDual : iterate)(D,n,Bindices); if(case != OPTIMAL) @@ -257,9 +284,8 @@ struct simplex { x[j]=0; for(int k=0; k < m; ++k) - x[Bindices[k]]=D[k][n]; - - cost=-Dm[n]; + 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 |