summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/base
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2019-03-27 22:27:57 +0000
committerKarl Berry <karl@freefriends.org>2019-03-27 22:27:57 +0000
commit9121faed0cf54a08882acf41120cf28fc3a32998 (patch)
treefe7f0c5e9ed02fe4a1fc389c27c88750c2807e43 /Build/source/utils/asymptote/base
parentb8e24fbf964b4030f68a8484a53d405869693266 (diff)
asy 2.48 sources
git-svn-id: svn://tug.org/texlive/trunk@50622 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/utils/asymptote/base')
-rw-r--r--Build/source/utils/asymptote/base/plain_prethree.asy15
-rw-r--r--Build/source/utils/asymptote/base/rational.asy73
-rw-r--r--Build/source/utils/asymptote/base/rationalSimplex.asy222
-rw-r--r--Build/source/utils/asymptote/base/shaders/fragment.glsl85
-rw-r--r--Build/source/utils/asymptote/base/shaders/vertex.glsl34
-rw-r--r--Build/source/utils/asymptote/base/simplex.asy142
-rw-r--r--Build/source/utils/asymptote/base/three.asy29
-rw-r--r--Build/source/utils/asymptote/base/three_light.asy27
-rw-r--r--Build/source/utils/asymptote/base/three_surface.asy38
9 files changed, 485 insertions, 180 deletions
diff --git a/Build/source/utils/asymptote/base/plain_prethree.asy b/Build/source/utils/asymptote/base/plain_prethree.asy
index cfc0c3c7411..0d3921e3f0d 100644
--- a/Build/source/utils/asymptote/base/plain_prethree.asy
+++ b/Build/source/utils/asymptote/base/plain_prethree.asy
@@ -165,7 +165,6 @@ struct light {
real[][] specular;
pen background=nullpen; // Background color of the 3D canvas.
real specularfactor;
- bool viewport; // Are the lights specified (and fixed) in the viewport frame?
triple[] position; // Only directional lights are currently implemented.
transform3 T=identity(4); // Transform to apply to normal vectors.
@@ -176,7 +175,7 @@ struct light {
pen[] ambient=array(diffuse.length,black),
pen[] specular=diffuse, pen background=nullpen,
real specularfactor=1,
- bool viewport=false, triple[] position) {
+ triple[] position) {
int n=diffuse.length;
assert(ambient.length == n && specular.length == n && position.length == n);
@@ -192,21 +191,18 @@ struct light {
this.position[i]=unit(position[i]);
}
this.specularfactor=specularfactor;
- this.viewport=viewport;
}
void operator init(pen diffuse=white, pen ambient=black, pen specular=diffuse,
- pen background=nullpen, real specularfactor=1,
- bool viewport=false...triple[] position) {
+ pen background=nullpen, real specularfactor=1 ...triple[] position) {
int n=position.length;
operator init(array(n,diffuse),array(n,ambient),array(n,specular),
- background,specularfactor,viewport,position);
+ background,specularfactor,position);
}
void operator init(pen diffuse=white, pen ambient=black, pen specular=diffuse,
- pen background=nullpen, bool viewport=false,
- real x, real y, real z) {
- operator init(diffuse,ambient,specular,background,viewport,(x,y,z));
+ pen background=nullpen, real x, real y, real z) {
+ operator init(diffuse,ambient,specular,background,(x,y,z));
}
void operator init(explicit light light) {
@@ -215,7 +211,6 @@ struct light {
specular=copy(light.specular);
background=light.background;
specularfactor=light.specularfactor;
- viewport=light.viewport;
position=copy(light.position);
}
diff --git a/Build/source/utils/asymptote/base/rational.asy b/Build/source/utils/asymptote/base/rational.asy
index a142063806d..d24da0d6645 100644
--- a/Build/source/utils/asymptote/base/rational.asy
+++ b/Build/source/utils/asymptote/base/rational.asy
@@ -133,6 +133,26 @@ bool operator ==(rational[][] r, rational[][] s)
return all(sequence(new bool(int i) {return r[i] == s[i];},r.length));
}
+bool[] operator <(rational[] r, rational s)
+{
+ return sequence(new bool(int i) {return r[i] < s;},r.length);
+}
+
+bool[] operator >(rational[] r, rational s)
+{
+ return sequence(new bool(int i) {return r[i] > s;},r.length);
+}
+
+bool[] operator <=(rational[] r, rational s)
+{
+ return sequence(new bool(int i) {return r[i] <= s;},r.length);
+}
+
+bool[] operator >=(rational[] r, rational s)
+{
+ return sequence(new bool(int i) {return r[i] >= s;},r.length);
+}
+
rational min(rational a, rational b)
{
return a <= b ? a : b;
@@ -143,35 +163,47 @@ rational max(rational a, rational b)
return a >= b ? a : b;
}
+string string(rational r)
+{
+ return r.q == 1 ? string(r.p) : string(r.p)+"/"+string(r.q);
+}
+
+string texstring(rational r)
+{
+ if(r.q == 1) return string(r.p);
+ string s;
+ if(r.p < 0) s="-";
+ return s+"\frac{"+string(abs(r.p))+"}{"+string(r.q)+"}";
+}
-void write(string s="", rational r, suffix suffix=endl) {
- if(r.q == 1)
- write(s+string(r.p),suffix);
- else
- write(s+string(r.p)+"/"+string(r.q),suffix);
+void write(file fout=stdout, string s="", rational r, suffix suffix=none)
+{
+ write(fout,s+string(r),suffix);
}
-void write(string s="", rational[] a, suffix suffix=endl) {
+void write(file fout=stdout, string s="", rational[] a, suffix suffix=none)
+{
if(s != "")
- write(s);
+ write(fout,s);
for(int i=0; i < a.length; ++i) {
- write(i,none);
- write(':\t',a[i]);
+ write(fout,i,none);
+ write(fout,':\t',a[i],endl);
}
- write(suffix);
+ write(fout,suffix);
}
-void write(string s="", rational[][] a, suffix suffix=endl) {
+void write(file fout=stdout, string s="", rational[][] a, suffix suffix=none)
+{
if(s != "")
- write(s);
+ write(fout,s);
for(int i=0; i < a.length; ++i) {
rational[] ai=a[i];
for(int j=0; j < ai.length; ++j) {
- write(ai[j],tab);
+ write(fout,ai[j],tab);
}
- write();
+ write(fout,endl);
}
- write(suffix);
+ write(fout,suffix);
}
bool rectangular(rational[][] m)
@@ -203,9 +235,18 @@ rational max(rational[] a)
rational abs(rational r)
{
- return rational(abs(r.p),r.q,false);
+ return rational(abs(r.p),r.q,false);
+}
+
+rational[] operator -(rational[] r)
+{
+ return sequence(new rational(int i) {return -r[i];},r.length);
}
+rational[][] rationalidentity(int n)
+{
+ return sequence(new rational[](int i) {return sequence(new rational(int j) {return j == i ? 1 : 0;},n);},n);
+}
/*
rational r=rational(1,3)+rational(1,4);
diff --git a/Build/source/utils/asymptote/base/rationalSimplex.asy b/Build/source/utils/asymptote/base/rationalSimplex.asy
index fd7128b8025..5175b87ba74 100644
--- a/Build/source/utils/asymptote/base/rationalSimplex.asy
+++ b/Build/source/utils/asymptote/base/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);
-*/
diff --git a/Build/source/utils/asymptote/base/shaders/fragment.glsl b/Build/source/utils/asymptote/base/shaders/fragment.glsl
new file mode 100644
index 00000000000..b6a85597451
--- /dev/null
+++ b/Build/source/utils/asymptote/base/shaders/fragment.glsl
@@ -0,0 +1,85 @@
+struct Material
+{
+ vec4 diffuse,ambient,emissive,specular;
+ float shininess;
+};
+
+struct Light
+{
+ vec4 direction;
+ vec4 diffuse,ambient,specular;
+};
+
+uniform int nlights;
+
+uniform Light lights[Nlights];
+
+uniform MaterialBuffer {
+ Material Materials[Nmaterials];
+};
+
+in vec3 Normal;
+
+#ifdef EXPLICIT_COLOR
+in vec4 Color;
+#endif
+flat in int materialIndex;
+
+out vec4 outColor;
+
+void main()
+{
+ vec4 Diffuse;
+ vec4 Ambient;
+ vec4 Emissive;
+ vec4 Specular;
+ float Shininess;
+
+#ifdef EXPLICIT_COLOR
+ if(materialIndex < 0) {
+ int index=-materialIndex-1;
+ Material m=Materials[index];
+ Diffuse=Color;
+ Ambient=Color;
+ Emissive=vec4(0.0,0.0,0.0,1.0);
+ Specular=m.specular;
+ Shininess=m.shininess;
+ } else {
+ Material m=Materials[materialIndex];
+ Diffuse=m.diffuse;
+ Ambient=m.ambient;
+ Emissive=m.emissive;
+ Specular=m.specular;
+ Shininess=m.shininess;
+ }
+#else
+ Material m=Materials[materialIndex];
+ Diffuse=m.diffuse;
+ Ambient=m.ambient;
+ Emissive=m.emissive;
+ Specular=m.specular;
+ Shininess=m.shininess;
+#endif
+ // Phong-Blinn model
+ if(nlights > 0) {
+ vec3 diffuse=vec3(0,0,0);
+ vec3 specular=vec3(0,0,0);
+ vec3 ambient=vec3(0,0,0);
+ vec3 Z=vec3(0,0,1);
+
+ for(int i=0; i < nlights; ++i) {
+ vec3 L=normalize(lights[i].direction.xyz);
+ diffuse += lights[i].diffuse.rgb*abs(dot(Normal,L));
+ ambient += lights[i].ambient.rgb;
+ specular += pow(abs(dot(Normal,normalize(L+Z))),Shininess)*
+ lights[i].specular.rgb;
+ }
+
+ vec3 color=diffuse*Diffuse.rgb+
+ ambient*Ambient.rgb+
+ specular*Specular.rgb+
+ Emissive.rgb;
+ outColor=vec4(color,Diffuse[3]);
+ } else
+ outColor=Diffuse;
+}
diff --git a/Build/source/utils/asymptote/base/shaders/vertex.glsl b/Build/source/utils/asymptote/base/shaders/vertex.glsl
new file mode 100644
index 00000000000..9e3feda6988
--- /dev/null
+++ b/Build/source/utils/asymptote/base/shaders/vertex.glsl
@@ -0,0 +1,34 @@
+in vec3 position;
+in vec3 normal;
+
+#ifdef EXPLICIT_COLOR
+in uint color;
+#endif
+
+in int material;
+
+uniform mat4 projViewMat;
+uniform mat4 viewMat;
+uniform mat4 normMat;
+
+out vec3 ViewPosition;
+out vec3 Normal;
+
+#ifdef EXPLICIT_COLOR
+out vec4 Color;
+#endif
+
+flat out int materialIndex;
+
+void main()
+{
+ gl_Position=projViewMat*vec4(position,1.0);
+ ViewPosition=(viewMat*vec4(position,1.0)).xyz;
+ Normal=normalize((normMat*vec4(normal,0)).xyz);
+
+#ifdef EXPLICIT_COLOR
+ Color=unpackUnorm4x8(color);
+#endif
+
+ materialIndex=material;
+}
diff --git a/Build/source/utils/asymptote/base/simplex.asy b/Build/source/utils/asymptote/base/simplex.asy
index 16a495a082e..040621dd904 100644
--- a/Build/source/utils/asymptote/base/simplex.asy
+++ b/Build/source/utils/asymptote/base/simplex.asy
@@ -1,4 +1,4 @@
-// General simplex solver written by John C. Bowman and Pouria Ramazi, 2018.
+// Real simplex solver written by John C. Bowman and Pouria Ramazi, 2018.
struct simplex {
static int OPTIMAL=0;
@@ -49,7 +49,7 @@ struct simplex {
if(Em[J] < 0) break;
if(J == N)
- return 0;
+ break;
int I=-1;
real M;
@@ -65,28 +65,69 @@ struct simplex {
real e=E[i][J];
if(e > epsilonA) {
real 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.
+ // Generate new tableau
Bindices[I]=J;
+ rowreduce(E,N,I,J);
+ }
+ return OPTIMAL;
+ }
+
+ int iterateDual(real[][] E, int N, int[] Bindices) {
+ while(true) {
+ // Find first negative entry in right (basic variable) column
+ real[] Em=E[m];
+ int I;
+ for(I=0; I < m; ++I) {
+ if(E[I][N] < 0) break;
+ }
+
+ if(I == m)
+ break;
+
+ int J=-1;
+ real M;
+ for(int j=0; j < N; ++j) {
+ real e=E[I][j];
+ if(e < epsilonA) {
+ M=-E[m][j]/e;
+ J=j;
+ break;
+ }
+ }
+ for(int j=J+1; j < N; ++j) {
+ real e=E[I][j];
+ 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)
+ return INFEASIBLE; // Can only happen in Phase 2.
// 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 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, real[] b, bool phase1=true) {
+ // 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,
+ bool dual=false) {
+ if(dual) phase1=false;
static real epsilon=sqrt(realEpsilon);
epsilonA=epsilon*norm(A);
- // Phase 1
+ // Phase 1
m=A.length;
if(m == 0) {case=INFEASIBLE; return;}
n=A[0].length;
@@ -102,7 +143,7 @@ struct simplex {
for(int i=0; i < m; ++i) {
real[] Ai=A[i];
real[] Ei=E[i];
- if(b[i] >= 0) {
+ if(b[i] >= 0 || dual) {
for(int j=0; j < n; ++j) {
real Aij=Ai[j];
Ei[j]=Aij;
@@ -130,7 +171,7 @@ struct simplex {
real sum=0;
for(int i=0; i < m; ++i) {
- real B=abs(b[i]);
+ real B=dual ? b[i] : abs(b[i]);
E[i][N]=B;
sum -= B;
}
@@ -140,27 +181,40 @@ struct simplex {
for(int j=0; j < m; ++j)
Em[n+j]=0.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(abs(Em[J]) > epsilonA) {
case=INFEASIBLE;
return;
}
- }
+ } else Bindices=sequence(new int(int x){return x;},m)+n-m;
+ real[] cB=phase1 ? new real[m] : c[n-m:n];
real[][] D=phase1 ? new real[m+1][n+1] : E;
- real[] Dm=D[m];
- real[] cb=phase1 ? new real[m] : c[n-m:n];
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=0; j < n; ++j)
+ if(Ei[j] != 0) 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];
+ cB[ip]=c[k];
real[] Dip=D[ip];
real[] Ei=E[i];
for(int j=0; j < n; ++j)
@@ -175,27 +229,29 @@ struct simplex {
Dip[j]=Em[j];
Dip[n]=Em[N];
- m=ip;
-
- for(int j=0; j < n; ++j) {
- real 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;
}
+ }
- // Done with Phase 1
+ real[] Dm=D[m];
+ for(int j=0; j < n; ++j) {
+ real sum=0;
+ for(int k=0; k < m; ++k)
+ sum += cB[k]*D[k][j];
+ Dm[j]=c[j]-sum;
}
-
+
real 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;
- }
+ case=(dual ? iterateDual : iterate)(D,n,Bindices);
+ if(case != OPTIMAL)
+ return;
for(int j=0; j < n; ++j)
x[j]=0;
@@ -204,7 +260,6 @@ 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
@@ -232,6 +287,9 @@ struct simplex {
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)
@@ -240,14 +298,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;
+ 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(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); // TODO: Check
- bool phase1=true;
- operator init(concat(c,array(count,0.0)),a,b,phase1);
+ operator init(concat(c,array(count,0.0)),a,b,phase1,dual);
- if(case == OPTIMAL)
+ if(case == OPTIMAL && count > 0)
x.delete(n,n+count-1);
}
}
diff --git a/Build/source/utils/asymptote/base/three.asy b/Build/source/utils/asymptote/base/three.asy
index 742ceda1c04..8b9d2137fc9 100644
--- a/Build/source/utils/asymptote/base/three.asy
+++ b/Build/source/utils/asymptote/base/three.asy
@@ -2021,7 +2021,7 @@ path3 arc(triple c, triple v1, triple v2, triple normal=O, bool direction=CCW)
}
string invalidnormal="invalid normal vector";
- real fuzz=sqrtEpsilon*max(abs(v1),abs(v2));
+ real fuzz=sqrtEpsilon;
if(abs(v1.z) > fuzz || abs(v2.z) > fuzz)
abort(invalidnormal);
@@ -2558,23 +2558,6 @@ private string Format(transform3 t, string sep=" ")
Format(t[0][3])+sep+Format(t[1][3])+sep+Format(t[2][3]);
}
-string lightscript(light light) {
- string script="for(var i=scene.lights.count-1; i >= 0; i--)
- scene.lights.removeByIndex(i);"+'\n\n';
- for(int i=0; i < light.position.length; ++i) {
- string Li="L"+string(i);
- real[] diffuse=light.diffuse[i];
- script += Li+"=scene.createLight();"+'\n'+
- Li+".direction.set("+format(-light.position[i],",")+");"+'\n'+
- Li+".color.set("+format((diffuse[0],diffuse[1],diffuse[2]),",")+");"+'\n';
- }
- // Work around initialization bug in Adobe Reader 8.0:
- return script +"
-scene.lightScheme=scene.LIGHT_MODE_HEADLAMP;
-scene.lightScheme=scene.LIGHT_MODE_FILE;
-";
-}
-
void writeJavaScript(string name, string preamble, string script)
{
file out=output(name);
@@ -2625,10 +2608,6 @@ string embed3D(string prefix, string label=prefix, string text=label,
shipout3(prefix,f);
string name=prefix+".js";
- // Adobe Reader doesn't appear to support user-specified viewport lights.
- bool lightscript=light.on() && !light.viewport;
- if(lightscript)
- writeJavaScript(name,lightscript(light),script);
if(!settings.inlinetex && !prconly())
file3.push(prefix+".prc");
@@ -2637,7 +2616,7 @@ string embed3D(string prefix, string label=prefix, string text=label,
transform3 inv=inverse(flipxz*P.T.modelview);
string options3="3Dlights="+
- (light.on() ? (light.viewport ? "Headlamp" : "File") : "None");
+ (light.on() ? "Headlamp" : "None");
if(defaultembed3Doptions != "") options3 += ","+defaultembed3Doptions;
if((settings.render < 0 || !settings.embed) && settings.auto3D)
@@ -2653,8 +2632,6 @@ string embed3D(string prefix, string label=prefix, string text=label,
if(options != "") options3 += ","+options;
if(settings.inlinetex)
prefix=jobname(prefix);
- if(lightscript)
- options3 += ",add3Djscript="+prefix+".js";
options3 += ",add3Djscript=asylabels.js";
return text == "" ? Embed(prefix+".prc","",options3,width,height) :
@@ -2909,7 +2886,7 @@ object embed(string prefix=outprefix(), string label=prefix,
P.zoom,m,M,P.viewportshift,
tinv*inv*shift(0,0,zcenter),Light.background(),Light.position,
Light.diffuse,Light.ambient,Light.specular,
- Light.viewport,view && !preview);
+ view && !preview);
if(!preview) return F;
}
diff --git a/Build/source/utils/asymptote/base/three_light.asy b/Build/source/utils/asymptote/base/three_light.asy
index 6cbd9b6f0be..3bd747dd4ec 100644
--- a/Build/source/utils/asymptote/base/three_light.asy
+++ b/Build/source/utils/asymptote/base/three_light.asy
@@ -86,36 +86,41 @@ pen color(triple normal, material m, light light, transform3 T=light.T) {
real[] Ambient=rgba(m.ambient());
real[] Specular=rgba(m.specular());
real[] p=rgba(m.emissive());
+ real[] diffuse={0,0,0,0};
+ real[] ambient={0,0,0,0};
+ real[] specular={0,0,0,0};
for(int i=0; i < position.length; ++i) {
- triple L=light.viewport ? position[i] : T*position[i];
- real Ldotn=max(dot(normal,L),0);
- p += light.ambient[i]*Ambient+Ldotn*light.diffuse[i]*Diffuse;
- // Apply specularfactor to partially compensate non-pixel-based rendering.
- if(Ldotn > 0) // Phong-Blinn model of specular reflection
- p += dot(normal,unit(L+Z))^s*light.specularfactor*
- light.specular[i]*Specular;
+ triple L=position[i];
+ real dotproduct=abs(dot(normal,L));
+ diffuse += dotproduct*light.diffuse[i];
+ ambient += light.ambient[i];
+ dotproduct=abs(dot(normal,unit(L+Z)));
+ // Phong-Blinn model of specular reflection
+ specular += dotproduct^s*light.specular[i];
}
+ p += diffuse*Diffuse;
+ p += ambient*Ambient;
+ // Apply specularfactor to partially compensate non-pixel-based rendering.
+ p += specular*Specular*light.specularfactor;
return rgb(p[0],p[1],p[2])+opacity(opacity(m.diffuse()));
}
light operator * (transform3 t, light light)
{
light light=light(light);
- if(!light.viewport) light.position=shiftless(t)*light.position;
return light;
}
light operator cast(triple v) {return light(v);}
-light Viewport=light(ambient=gray(0.1),specularfactor=3,viewport=true,
- (0.25,-0.25,1));
+light Viewport=light(ambient=gray(0.1),specularfactor=3,(0.25,-0.25,1));
light White=light(new pen[] {rgb(0.38,0.38,0.45),rgb(0.6,0.6,0.67),
rgb(0.5,0.5,0.57)},specularfactor=3,
new triple[] {(-2,-1.5,-0.5),(2,1.1,-2.5),(-0.5,0,2)});
light Headlamp=light(gray(0.8),ambient=gray(0.1),specular=gray(0.7),
- specularfactor=3,viewport=true,dir(42,48));
+ specularfactor=3,dir(42,48));
currentlight=Headlamp;
diff --git a/Build/source/utils/asymptote/base/three_surface.asy b/Build/source/utils/asymptote/base/three_surface.asy
index 6e10246c774..b3eebee45cc 100644
--- a/Build/source/utils/asymptote/base/three_surface.asy
+++ b/Build/source/utils/asymptote/base/three_surface.asy
@@ -1409,11 +1409,32 @@ void draw3D(frame f, int type=0, patch s, triple center=O, material m,
PRCshininess,s.colors,interaction.type,prc);
}
+int computeNormals(triple[] v, int[][] vi, triple[] n, int[][] ni)
+{
+ triple lastnormal=O;
+ for(int i=0; i < vi.length; ++i) {
+ int[] vii=vi[i];
+ int[] nii=ni[i];
+ triple normal=normal(new triple[] {v[vii[0]],v[vii[1]],v[vii[2]]});
+ if(normal != lastnormal || n.length == 0) {
+ n.push(normal);
+ lastnormal=normal;
+ }
+ nii[0]=nii[1]=nii[2]=n.length-1;
+ }
+ return ni.length;
+}
+
// Draw triangles on a frame.
void draw(frame f, triple[] v, int[][] vi,
triple[] n={}, int[][] ni={}, material m=currentpen, pen[] p={},
int[][] pi={}, light light=currentlight)
{
+ bool normals=ni.length > 0;
+ if(!normals) {
+ ni=new int[vi.length][3];
+ normals=computeNormals(v,vi,n,ni) > 0;
+ }
if(p.length > 0)
m=mean(p);
m=material(m,light);
@@ -1428,23 +1449,12 @@ void draw(picture pic=currentpicture, triple[] v, int[][] vi,
triple[] n={}, int[][] ni={}, material m=currentpen, pen[] p={},
int[][] pi={}, light light=currentlight)
{
- bool colors=pi.length > 0;
bool normals=ni.length > 0;
- if(!colors && !normals) {
- n=new triple[];
+ if(!normals) {
ni=new int[vi.length][3];
- triple lastnormal=O;
- for(int i=0; i < vi.length; ++i) {
- int[] vii=vi[i];
- int[] nii=ni[i];
- triple normal=normal(new triple[] {v[vii[0]],v[vii[1]],v[vii[2]]});
- if(normal != lastnormal || n.length == 0) {
- n.push(normal);
- lastnormal=normal;
- }
- nii[0]=nii[1]=nii[2]=n.length-1;
- }
+ normals=computeNormals(v,vi,n,ni) > 0;
}
+ bool colors=pi.length > 0;
pic.add(new void(frame f, transform3 t, picture pic, projection P) {
triple[] v=t*v;