summaryrefslogtreecommitdiff
path: root/graphics/asymptote/base/simplex.asy
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2019-11-07 03:01:39 +0000
committerNorbert Preining <norbert@preining.info>2019-11-07 03:01:39 +0000
commit2e1d63b8ed8c6b7c6d206bfe9e2712797108e8bd (patch)
tree7ea198b043998590117b2730e9dec10327ef3da3 /graphics/asymptote/base/simplex.asy
parent590fd8b560523cdaea19c59aea61e781138e87f8 (diff)
CTAN sync 201911070301
Diffstat (limited to 'graphics/asymptote/base/simplex.asy')
-rw-r--r--graphics/asymptote/base/simplex.asy156
1 files changed, 91 insertions, 65 deletions
diff --git a/graphics/asymptote/base/simplex.asy b/graphics/asymptote/base/simplex.asy
index 040621dd90..b9ff172077 100644
--- a/graphics/asymptote/base/simplex.asy
+++ b/graphics/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