summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/asymptote/rationalSimplex.asy
diff options
context:
space:
mode:
Diffstat (limited to 'Master/texmf-dist/asymptote/rationalSimplex.asy')
-rw-r--r--Master/texmf-dist/asymptote/rationalSimplex.asy222
1 files changed, 152 insertions, 70 deletions
diff --git a/Master/texmf-dist/asymptote/rationalSimplex.asy b/Master/texmf-dist/asymptote/rationalSimplex.asy
index fd7128b8025..5175b87ba74 100644
--- a/Master/texmf-dist/asymptote/rationalSimplex.asy
+++ b/Master/texmf-dist/asymptote/rationalSimplex.asy
@@ -1,6 +1,29 @@
// Rational simplex solver written by John C. Bowman and Pouria Ramazi, 2018.
import rational;
+void simplexTableau(rational[][] E, int[] Bindices, int I=-1, int J=-1) {}
+void simplexPhase2() {}
+
+void simplexWrite(rational[][] E, int[] Bindicies, int, int)
+{
+ int m=E.length-1;
+ int n=E[0].length-1;
+
+ write(E[m][n],tab);
+ for(int j=0; j < n; ++j)
+ write(E[m][j],tab);
+ write();
+
+ for(int i=0; i < m; ++i) {
+ write(E[i][n],tab);
+ for(int j=0; j < n; ++j) {
+ write(E[i][j],tab);
+ }
+ write();
+ }
+ write();
+};
+
struct simplex {
static int OPTIMAL=0;
static int UNBOUNDED=1;
@@ -49,7 +72,7 @@ struct simplex {
if(Em[J] < 0) break;
if(J == N)
- return 0;
+ break;
int I=-1;
rational M;
@@ -65,28 +88,74 @@ struct simplex {
rational e=E[i][J];
if(e > 0) {
rational v=E[i][N]/e;
- if(v <= M) {M=v; I=i;}
+ if(v < M) {M=v; I=i;} // Bland's rule: choose smallest argmin
}
}
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);
+ }
+ return OPTIMAL;
+ }
+
+ int iterateDual(rational[][] E, int N, int[] Bindices) {
+ while(true) {
+ // Find first negative entry in right (basic variable) column
+ rational[] Em=E[m];
+ int I;
+ for(I=0; I < m; ++I) {
+ if(E[I][N] < 0) break;
+ }
+
+ if(I == m)
+ break;
+
+ int J=-1;
+ rational M;
+ for(int j=0; j < N; ++j) {
+ rational e=E[I][j];
+ if(e < 0) {
+ M=-E[m][j]/e;
+ J=j;
+ break;
+ }
+ }
+ for(int j=J+1; j < N; ++j) {
+ rational e=E[I][j];
+ if(e < 0) {
+ rational v=-E[m][j]/e;
+ if(v < M) {M=v; J=j;} // Bland's rule: choose smallest argmin
+ }
+ }
+ if(J == -1)
+ 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 0;
+ 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 length n, b is a
- // vector of length m, and c is a vector of length n.
+ // 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
+ bool phase1=true, bool dual=false) {
+ if(dual) phase1=false;
+ // Phase 1
m=A.length;
+ if(m == 0) {case=INFEASIBLE; return;}
n=A[0].length;
+ if(n == 0) {case=INFEASIBLE; return;}
int N=phase1 ? n+m : n;
rational[][] E=new rational[m+1][N+1];
@@ -98,7 +167,7 @@ struct simplex {
for(int i=0; i < m; ++i) {
rational[] Ai=A[i];
rational[] Ei=E[i];
- if(b[i] >= 0) {
+ if(b[i] >= 0 || dual) {
for(int j=0; j < n; ++j) {
rational Aij=Ai[j];
Ei[j]=Aij;
@@ -126,7 +195,7 @@ struct simplex {
rational sum=0;
for(int i=0; i < m; ++i) {
- rational B=abs(b[i]);
+ rational B=dual ? b[i] : abs(b[i]);
E[i][N]=B;
sum -= B;
}
@@ -136,27 +205,45 @@ struct simplex {
for(int j=0; j < m; ++j)
Em[n+j]=0;
- int[] Bindices=sequence(new int(int x){return x;},m)+n;
+ int[] Bindices;
if(phase1) {
+ Bindices=sequence(new int(int x){return x;},m)+n;
iterate(E,N,Bindices);
if(Em[J] != 0) {
+ simplexTableau(E,Bindices);
case=INFEASIBLE;
return;
}
- }
-
+ } else Bindices=sequence(new int(int x){return x;},m)+n-m;
+
+ rational[] cB=phase1 ? new rational[m] : c[n-m:n];
rational[][] D=phase1 ? new rational[m+1][n+1] : E;
- rational[] Dm=D[m];
- rational[] cb=phase1 ? new rational[m] : c[n-m:n];
if(phase1) {
+ bool output=true;
+ // Drive artificial variables out of basis.
+ for(int i=0; i < m; ++i) {
+ int k=Bindices[i];
+ if(k >= n) {
+ rational[] Ei=E[i];
+ int j;
+ for(j=0; j < n; ++j)
+ if(Ei[j] != 0) break;
+ if(j == n) continue;
+ output=false;
+ simplexTableau(E,Bindices,i,j);
+ Bindices[i]=j;
+ rowreduce(E,n,i,j);
+ }
+ }
+ if(output) simplexTableau(E,Bindices);
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];
+ cB[ip]=c[k];
rational[] Dip=D[ip];
rational[] Ei=E[i];
for(int j=0; j < n; ++j)
@@ -171,27 +258,33 @@ struct simplex {
Dip[j]=Em[j];
Dip[n]=Em[N];
- m=ip;
-
- for(int j=0; j < n; ++j) {
- rational sum=0;
- for(int k=0; k < m; ++k)
- sum += cb[k]*D[k][j];
- Dm[j]=c[j]-sum;
+ if(m > ip) {
+ Bindices.delete(ip,m-1);
+ D.delete(ip,m-1);
+ m=ip;
}
+ if(!output) simplexTableau(D,Bindices);
+ }
- // Done with Phase 1
+ rational[] Dm=D[m];
+ for(int j=0; j < n; ++j) {
+ rational sum=0;
+ for(int k=0; k < m; ++k)
+ sum += cB[k]*D[k][j];
+ Dm[j]=c[j]-sum;
}
-
+
rational sum=0;
for(int k=0; k < m; ++k)
- sum += cb[k]*D[k][n];
+ sum += cB[k]*D[k][n];
Dm[n]=-sum;
- if(iterate(D,n,Bindices) == UNBOUNDED) {
- case=UNBOUNDED;
- return;
- }
+ simplexPhase2();
+
+ case=(dual ? iterateDual : iterate)(D,n,Bindices);
+ simplexTableau(D,Bindices);
+ if(case != OPTIMAL)
+ return;
for(int j=0; j < n; ++j)
x[j]=0;
@@ -200,15 +293,16 @@ struct simplex {
x[Bindices[k]]=D[k][n];
cost=-Dm[n];
- case=OPTIMAL;
}
// 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 length n, b is a
- // vector of length m, and c is a vector of length n.
+ // 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)
@@ -226,6 +320,9 @@ struct simplex {
int k=0;
+ bool phase1=false;
+ bool dual=count == m && all(c >= 0);
+
for(int i=0; i < m; ++i) {
rational[] ai=a[i];
for(int j=0; j < k; ++j)
@@ -234,47 +331,32 @@ struct simplex {
ai[n+k]=-s[i];
for(int j=k+1; j < count; ++j)
ai[n+j]=0;
- if(s[i] != 0) ++k;
+ int si=s[i];
+ 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(si*bi > 0) {
+ if(dual && si == 1) {
+ b[i]=-bi;
+ s[i]=-1;
+ for(int j=0; j < n+count; ++j)
+ ai[j]=-ai[j];
+ } else
+ phase1=true;
+ }
+ }
}
- bool phase1=!all(s == -1);
- operator init(concat(c,array(count,rational(0))),a,b,phase1);
+ operator init(concat(c,array(count,rational(0))),a,b,phase1,dual);
- if(case == OPTIMAL)
+ if(case == OPTIMAL && count > 0)
x.delete(n,n+count-1);
}
}
-
-/*
-simplex S=simplex(new rational[] {4,1,1},
- new rational[][] {{2,1,2},{3,3,1}},
- new rational[] {4,3});
-
-simplex S=simplex(new rational[] {2,6,1,1},
- new rational[][] {{1,2,0,1},{1,2,1,1},{1,3,-1,2},{1,1,1,0}},
- new rational[] {6,7,7,5});
-simplex S=simplex(new rational[] {-10,-12,-12,0,0,0},
- new rational[][] {{1,2,2,1,0,0},
- {2,1,2,0,1,0},
- {2,2,1,0,0,1}},
- new rational[] {20,20,20});
-
-simplex S=simplex(new rational[] {-10,-12,-12},
- new rational[][] {{1,2,2},
- {2,1,2},
- {2,2,1}},
- new int[] {0,0,-1},
- new rational[] {20,20,20});
-
-simplex S=simplex(new rational[] {1,1,1,0},
- new rational[][] {{1,2,3,0},
- {-1,2,6,0},
- {0,4,9,0},
- {0,0,3,1}},
- new rational[] {3,2,5,1});
-
-write();
-write("case:",S.case);
-write("x:",S.x);
-write("Cost=",S.cost);
-*/