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.asy280
1 files changed, 280 insertions, 0 deletions
diff --git a/Master/texmf-dist/asymptote/rationalSimplex.asy b/Master/texmf-dist/asymptote/rationalSimplex.asy
new file mode 100644
index 00000000000..fd7128b8025
--- /dev/null
+++ b/Master/texmf-dist/asymptote/rationalSimplex.asy
@@ -0,0 +1,280 @@
+// Rational simplex solver written by John C. Bowman and Pouria Ramazi, 2018.
+import rational;
+
+struct simplex {
+ static int OPTIMAL=0;
+ static int UNBOUNDED=1;
+ static int INFEASIBLE=2;
+
+ int case;
+ rational[] x;
+ rational cost;
+
+ int m,n;
+ int J;
+
+ // Row reduce based on pivot E[I][J]
+ void rowreduce(rational[][] E, int N, int I, int J) {
+ rational[] EI=E[I];
+ rational v=EI[J];
+ for(int j=0; j < J; ++j) EI[j] /= v;
+ EI[J]=1;
+ for(int j=J+1; j <= N; ++j) EI[j] /= v;
+
+ for(int i=0; i < I; ++i) {
+ rational[] Ei=E[i];
+ rational EiJ=Ei[J];
+ for(int j=0; j < J; ++j)
+ Ei[j] -= EI[j]*EiJ;
+ Ei[J]=0;
+ for(int j=J+1; j <= N; ++j)
+ Ei[j] -= EI[j]*EiJ;
+ }
+ for(int i=I+1; i <= m; ++i) {
+ rational[] Ei=E[i];
+ rational EiJ=Ei[J];
+ for(int j=0; j < J; ++j)
+ Ei[j] -= EI[j]*EiJ;
+ Ei[J]=0;
+ for(int j=J+1; j <= N; ++j)
+ Ei[j] -= EI[j]*EiJ;
+ }
+ }
+
+ int iterate(rational[][] E, int N, int[] Bindices) {
+ while(true) {
+ // Find first negative entry in bottom (reduced cost) row
+ rational[] Em=E[m];
+ for(J=0; J < N; ++J)
+ if(Em[J] < 0) break;
+
+ if(J == N)
+ return 0;
+
+ int I=-1;
+ rational M;
+ for(int i=0; i < m; ++i) {
+ rational e=E[i][J];
+ if(e > 0) {
+ M=E[i][N]/e;
+ I=i;
+ break;
+ }
+ }
+ for(int i=I+1; i < m; ++i) {
+ rational e=E[i][J];
+ if(e > 0) {
+ rational v=E[i][N]/e;
+ if(v <= M) {M=v; I=i;}
+ }
+ }
+ if(I == -1)
+ return UNBOUNDED; // Can only happen in Phase 2.
+
+ Bindices[I]=J;
+
+ // Generate new tableau
+ rowreduce(E,N,I,J);
+ }
+ return 0;
+ }
+
+ // 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.
+ void operator init(rational[] c, rational[][] A, rational[] b,
+ bool phase1=true) {
+ // Phase 1
+ m=A.length;
+ n=A[0].length;
+
+ int N=phase1 ? n+m : n;
+ rational[][] E=new rational[m+1][N+1];
+ rational[] Em=E[m];
+
+ for(int j=0; j < n; ++j)
+ Em[j]=0;
+
+ for(int i=0; i < m; ++i) {
+ rational[] Ai=A[i];
+ rational[] Ei=E[i];
+ if(b[i] >= 0) {
+ for(int j=0; j < n; ++j) {
+ rational Aij=Ai[j];
+ Ei[j]=Aij;
+ Em[j] -= Aij;
+ }
+ } else {
+ for(int j=0; j < n; ++j) {
+ rational Aij=-Ai[j];
+ Ei[j]=Aij;
+ Em[j] -= Aij;
+ }
+ }
+ }
+
+ if(phase1) {
+ for(int i=0; i < m; ++i) {
+ rational[] Ei=E[i];
+ for(int j=0; j < i; ++j)
+ Ei[n+j]=0;
+ Ei[n+i]=1;
+ for(int j=i+1; j < m; ++j)
+ Ei[n+j]=0;
+ }
+ }
+
+ rational sum=0;
+ for(int i=0; i < m; ++i) {
+ rational B=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;
+
+ int[] Bindices=sequence(new int(int x){return x;},m)+n;
+
+ if(phase1) {
+ iterate(E,N,Bindices);
+
+ if(Em[J] != 0) {
+ case=INFEASIBLE;
+ return;
+ }
+ }
+
+ 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) {
+ 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];
+ rational[] Dip=D[ip];
+ rational[] Ei=E[i];
+ for(int j=0; j < n; ++j)
+ Dip[j]=Ei[j];
+ Dip[n]=Ei[N];
+ ++ip;
+ }
+
+ rational[] Dip=D[ip];
+ rational[] Em=E[m];
+ for(int j=0; j < n; ++j)
+ 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;
+ }
+
+ // Done with Phase 1
+ }
+
+ rational sum=0;
+ for(int k=0; k < m; ++k)
+ sum += cb[k]*D[k][n];
+ Dm[n]=-sum;
+
+ if(iterate(D,n,Bindices) == UNBOUNDED) {
+ case=UNBOUNDED;
+ return;
+ }
+
+ for(int j=0; j < n; ++j)
+ x[j]=0;
+
+ for(int k=0; k < m; ++k)
+ 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.
+ void operator init(rational[] c, rational[][] A, int[] s, rational[] b) {
+ int m=A.length;
+ int n=A[0].length;
+
+ int count=0;
+ for(int i=0; i < m; ++i)
+ if(s[i] != 0) ++count;
+
+ rational[][] a=new rational[m][n+count];
+
+ for(int i=0; i < m; ++i) {
+ rational[] ai=a[i];
+ rational[] Ai=A[i];
+ for(int j=0; j < n; ++j) {
+ ai[j]=Ai[j];
+ }
+ }
+
+ int k=0;
+
+ for(int i=0; i < m; ++i) {
+ rational[] ai=a[i];
+ for(int j=0; j < k; ++j)
+ ai[n+j]=0;
+ if(k < count)
+ ai[n+k]=-s[i];
+ for(int j=k+1; j < count; ++j)
+ ai[n+j]=0;
+ if(s[i] != 0) ++k;
+ }
+
+ bool phase1=!all(s == -1);
+ operator init(concat(c,array(count,rational(0))),a,b,phase1);
+
+ if(case == OPTIMAL)
+ 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);
+*/