From 8d656f7b22badc7a1e4e48811521048ef3df7101 Mon Sep 17 00:00:00 2001 From: Karl Berry Date: Thu, 10 Jan 2019 18:49:48 +0000 Subject: asy 2.47 sources git-svn-id: svn://tug.org/texlive/trunk@49658 c570f23f-e606-0410-a88d-b1316a301751 --- Build/source/utils/asymptote/base/animation.asy | 4 - Build/source/utils/asymptote/base/asymptote.py | 10 +- Build/source/utils/asymptote/base/bsp.asy | 26 +- Build/source/utils/asymptote/base/plain.asy | 12 +- Build/source/utils/asymptote/base/plain_Label.asy | 13 +- Build/source/utils/asymptote/base/plain_boxes.asy | 4 +- .../utils/asymptote/base/plain_constants.asy | 3 + Build/source/utils/asymptote/base/plain_pens.asy | 8 +- .../source/utils/asymptote/base/plain_picture.asy | 92 ++-- .../source/utils/asymptote/base/plain_scaling.asy | 36 +- .../source/utils/asymptote/base/plain_shipout.asy | 49 ++- Build/source/utils/asymptote/base/plain_xasy.asy | 132 ------ Build/source/utils/asymptote/base/rational.asy | 5 + .../utils/asymptote/base/rationalSimplex.asy | 280 +++++++++++++ Build/source/utils/asymptote/base/simplex.asy | 463 ++++++++++----------- Build/source/utils/asymptote/base/slide.asy | 2 +- Build/source/utils/asymptote/base/three.asy | 7 +- 17 files changed, 682 insertions(+), 464 deletions(-) delete mode 100644 Build/source/utils/asymptote/base/plain_xasy.asy create mode 100644 Build/source/utils/asymptote/base/rationalSimplex.asy (limited to 'Build/source/utils/asymptote/base') diff --git a/Build/source/utils/asymptote/base/animation.asy b/Build/source/utils/asymptote/base/animation.asy index 270adce713a..e2fd1ea18c9 100644 --- a/Build/source/utils/asymptote/base/animation.asy +++ b/Build/source/utils/asymptote/base/animation.asy @@ -56,7 +56,6 @@ struct animation { string format=nativeformat(); plain.shipout(name,f,format=format,view=false); files.push(name+"."+format); - shipped=false; } void add(picture pic=currentpicture, enclosure enclosure=NoBox) { @@ -119,7 +118,6 @@ struct animation { plain.shipout(prefix,multi,view=view); settings.inlinetex=inlinetex; } - shipped=true; } string load(int frames, real delay=animationdelay, string options="", @@ -152,8 +150,6 @@ struct animation { if(global) export(filename,enclosure,multipage=multipage); - shipped=false; - if(!keep) { exitfcn currentexitfunction=atexit(); void exitfunction() { diff --git a/Build/source/utils/asymptote/base/asymptote.py b/Build/source/utils/asymptote/base/asymptote.py index 71e7b9beed6..19b72c50826 100755 --- a/Build/source/utils/asymptote/base/asymptote.py +++ b/Build/source/utils/asymptote/base/asymptote.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python3 + # Python module to feed Asymptote with commands # (modified from gnuplot.py) from subprocess import * @@ -6,7 +8,7 @@ class asy: self.session = Popen(['asy','-quiet','-inpipe=0','-outpipe=2'],stdin=PIPE) self.help() def send(self, cmd): - self.session.stdin.write(cmd+'\n') + self.session.stdin.write(bytes(cmd+'\n','utf-8')) self.session.stdin.flush() def size(self, size): self.send("size(%d);" % size) @@ -23,10 +25,10 @@ class asy: def erase(self): self.send("erase();") def help(self): - print "Asymptote session is open. Available methods are:" - print " help(), size(int), draw(str), fill(str), clip(str), label(str), shipout(str), send(str), erase()" + print("Asymptote session is open. Available methods are:") + print(" help(), size(int), draw(str), fill(str), clip(str), label(str), shipout(str), send(str), erase()") def __del__(self): - print "closing Asymptote session..." + print("closing Asymptote session...") self.send('quit'); self.session.stdin.close(); self.session.wait() diff --git a/Build/source/utils/asymptote/base/bsp.asy b/Build/source/utils/asymptote/base/bsp.asy index 2a4d93daa81..cf9d7e0781f 100644 --- a/Build/source/utils/asymptote/base/bsp.asy +++ b/Build/source/utils/asymptote/base/bsp.asy @@ -184,20 +184,20 @@ void add(picture pic=currentpicture, face[] faces, for(int i=0; i < n; ++i) Faces[i]=faces[i].copy(); - pic.nodes.push(new void (frame f, transform t, transform T, - pair m, pair M) { - // Fit all of the pictures so we know their exact sizes. - face[] faces=new face[n]; - for(int i=0; i < n; ++i) { - faces[i]=Faces[i].copy(); - face F=faces[i]; - F.t=t*T*F.pic.T; - F.fit=F.pic.fit(t,T*F.pic.T,m,M); - } + pic.add(new void (frame f, transform t, transform T, + pair m, pair M) { + // Fit all of the pictures so we know their exact sizes. + face[] faces=new face[n]; + for(int i=0; i < n; ++i) { + faces[i]=Faces[i].copy(); + face F=faces[i]; + F.t=t*T*F.pic.T; + F.fit=F.pic.fit(t,T*F.pic.T,m,M); + } - bsp bsp=bsp(faces,P); - if(bsp != null) bsp.add(f); - }); + bsp bsp=bsp(faces,P); + if(bsp != null) bsp.add(f); + }); for(int i=0; i < n; ++i) { picture F=Faces[i].pic; diff --git a/Build/source/utils/asymptote/base/plain.asy b/Build/source/utils/asymptote/base/plain.asy index 3bb8cf7ce93..4afe3cb3de4 100644 --- a/Build/source/utils/asymptote/base/plain.asy +++ b/Build/source/utils/asymptote/base/plain.asy @@ -41,18 +41,19 @@ include plain_debugger; typedef void exitfcn(); -bool needshipout() { - return !shipped && !currentpicture.empty(); -} - void updatefunction() { + implicitshipout=true; if(!currentpicture.uptodate) shipout(); + implicitshipout=false; } void exitfunction() { - if(needshipout()) shipout(); + implicitshipout=true; + if(!currentpicture.empty()) + shipout(); + implicitshipout=false; } atupdate(updatefunction); @@ -273,7 +274,6 @@ if(settings.autoimport != "") { string s=settings.autoimport; settings.autoimport=""; eval("import \""+s+"\" as dummy",true); - shipped=false; atupdate(updatefunction); atexit(exitfunction); settings.autoimport=s; diff --git a/Build/source/utils/asymptote/base/plain_Label.asy b/Build/source/utils/asymptote/base/plain_Label.asy index 516ed4be711..c50ce003be0 100644 --- a/Build/source/utils/asymptote/base/plain_Label.asy +++ b/Build/source/utils/asymptote/base/plain_Label.asy @@ -339,8 +339,8 @@ struct Label { } pic.add(new void (frame f, transform t) { - out(f,t,point(g,position), - alignrelative ? -Align*dir(t*g,position)*I : Align); + out(f,t,point(g,position),alignrelative ? + inverse(rotation(t))*-Align*dir(t*g,position)*I : Align); },!alignrelative); frame f; @@ -485,6 +485,15 @@ void label(picture pic=currentpicture, Label L, align align=NoAlign, { label(pic,L,L.position,align,p,filltype); } + +// Label, but with postscript coords instead of asy +void label(pair origin, picture pic=currentpicture, Label L, align align=NoAlign, + pen p=currentpen, filltype filltype=NoFill) +{ + picture opic; + label(opic,L,L.position,align,p,filltype); + add(pic,opic,origin); +} void label(picture pic=currentpicture, Label L, explicit path g, align align=NoAlign, pen p=currentpen, filltype filltype=NoFill) diff --git a/Build/source/utils/asymptote/base/plain_boxes.asy b/Build/source/utils/asymptote/base/plain_boxes.asy index 0c6401311d8..c675ab4cde9 100644 --- a/Build/source/utils/asymptote/base/plain_boxes.asy +++ b/Build/source/utils/asymptote/base/plain_boxes.asy @@ -131,7 +131,9 @@ frame bbox(picture pic=currentpicture, real xmargin=0, real ymargin=xmargin, pen p=currentpen, filltype filltype=NoFill) { - frame f=pic.fit(max(pic.xsize-2*xmargin,0),max(pic.ysize-2*ymargin,0)); + real penwidth=linewidth(p); + frame f=pic.fit(max(pic.xsize-2*(xmargin+penwidth),0), + max(pic.ysize-2*(ymargin+penwidth),0)); box(f,xmargin,ymargin,p,filltype,above=false); return f; } diff --git a/Build/source/utils/asymptote/base/plain_constants.asy b/Build/source/utils/asymptote/base/plain_constants.asy index bf6988d4c42..136b3540218 100644 --- a/Build/source/utils/asymptote/base/plain_constants.asy +++ b/Build/source/utils/asymptote/base/plain_constants.asy @@ -36,6 +36,9 @@ restricted real sqrtEpsilon=sqrt(realEpsilon); restricted pair Align=sqrtEpsilon*NE; restricted int mantissaBits=ceil(-log(realEpsilon)/log(2))+1; +restricted transform identity; +restricted transform zeroTransform=(0,0,0,0,0,0); + int min(... int[] a) {return min(a);} int max(... int[] a) {return max(a);} diff --git a/Build/source/utils/asymptote/base/plain_pens.asy b/Build/source/utils/asymptote/base/plain_pens.asy index 627ce2e8c84..22cf5ca4817 100644 --- a/Build/source/utils/asymptote/base/plain_pens.asy +++ b/Build/source/utils/asymptote/base/plain_pens.asy @@ -324,10 +324,16 @@ pen rgba(real[] a) // Return a pen corresponding to a given 6-character RGB hexidecimal string. pen rgb(string s) { - real value(string s, int i) {return hex(substr(s,2i,2))/255;} + int offset=substr(s,0,1) == '#' ? 1 : 0; + real value(string s, int i) {return hex(substr(s,2i+offset,2))/255;} return rgb(value(s,0),value(s,1),value(s,2)); } +pen RGB(int r, int g, int b) +{ + return rgb(r/255,g/255,b/255); +} + pen[] operator +(pen[] a, pen b) { return sequence(new pen(int i) {return a[i]+b;},a.length); diff --git a/Build/source/utils/asymptote/base/plain_picture.asy b/Build/source/utils/asymptote/base/plain_picture.asy index f1dccbb8bea..025e6741cfd 100644 --- a/Build/source/utils/asymptote/base/plain_picture.asy +++ b/Build/source/utils/asymptote/base/plain_picture.asy @@ -169,6 +169,15 @@ typedef void drawer(frame f, transform t); // TODO: Add documentation as to what T is. typedef void drawerBound(frame f, transform t, transform T, pair lb, pair rt); +struct node { + drawerBound d; + string key; + void operator init(drawerBound d, string key=xasyKEY()) { + this.d=d; + this.key=key; + } +} + // PairOrTriple <<<1 // This struct is used to represent a userMin/userMax which serves as both a // pair and a triple depending on the context. @@ -199,9 +208,18 @@ struct picture { // <<<1 typedef void drawerBound3(frame f, transform3 t, transform3 T, picture pic, projection P, triple lb, triple rt); + struct node3 { + drawerBound3 d; + string key; + void operator init(drawerBound3 d, string key=xasyKEY()) { + this.d=d; + this.key=key; + } + } + // The functions to do the deferred drawing. - drawerBound[] nodes; - drawerBound3[] nodes3; + node[] nodes; + node3[] nodes3; bool uptodate=true; @@ -266,15 +284,20 @@ struct picture { // <<<1 } init(); - // Erase the current picture, retaining any size specification. - void erase() { + // Erase the current picture, retaining bounds. + void clear() { nodes.delete(); nodes3.delete(); + legend.delete(); + } + + // Erase the current picture, retaining any size specification. + void erase() { + clear(); bounds.erase(); bounds3.erase(); T=identity(); scale=new ScaleT; - legend.delete(); init(); } @@ -445,15 +468,15 @@ struct picture { // <<<1 uptodate=false; if(!exact) bounds.exact=false; if(above) - nodes.push(d); + nodes.push(node(d)); else - nodes.insert(0,d); + nodes.insert(0,node(d)); } // Faster implementation of most common case. void addExactAbove(drawerBound d) { uptodate=false; - nodes.push(d); + nodes.push(node(d)); } void add(drawer d, bool exact=false, bool above=true) { @@ -466,9 +489,9 @@ struct picture { // <<<1 uptodate=false; if(!exact) bounds.exact=false; if(above) - nodes3.push(d); + nodes3.push(node3(d)); else - nodes3.insert(0,d); + nodes3.insert(0,node3(d)); } void add(drawer3 d, bool exact=false, bool above=true) { @@ -675,17 +698,20 @@ struct picture { // <<<1 frame fit(transform t, transform T0=T, pair m, pair M) { frame f; - int n = nodes.length; - for(int i=0; i < n; ++i) - nodes[i](f,t,T0,m,M); + for(node n : nodes) { + xasyKEY(n.key); + n.d(f,t,T0,m,M); + } return f; } frame fit3(transform3 t, transform3 T0=T3, picture pic, projection P, triple m, triple M) { frame f; - for(int i=0; i < nodes3.length; ++i) - nodes3[i](f,t,T0,pic,P,m,M); + for(node3 n : nodes3) { + xasyKEY(n.key); + n.d(f,t,T0,pic,P,m,M); + } return f; } @@ -905,7 +931,7 @@ struct picture { // <<<1 // picture. Fitting this picture will not scale as the original picture would. picture drawcopy() { picture dest=new picture; - dest.nodes = copy(nodes); + dest.nodes=copy(nodes); dest.nodes3=copy(nodes3); dest.T=T; dest.T3=T3; @@ -945,13 +971,15 @@ struct picture { // <<<1 picture dest=drawcopy(); // Replace nodes with a single drawer that realizes the transform. - drawerBound[] oldnodes = dest.nodes; + node[] oldnodes = dest.nodes; void drawAll(frame f, transform tt, transform T, pair lb, pair rt) { transform Tt = T*t; - for (var node : oldnodes) - node(f, tt, Tt, lb, rt); + for (node n : oldnodes) { + xasyKEY(n.key); + n.d(f,tt,Tt,lb,rt); + } } - dest.nodes = new drawerBound[] { drawAll }; + dest.nodes = new node[] {node(drawAll)}; dest.uptodate=uptodate; dest.bounds=bounds.transformed(t); @@ -984,16 +1012,19 @@ struct picture { // <<<1 picture srcCopy=src.drawcopy(); // Draw by drawing the copied picture. - if(srcCopy.nodes.length > 0) - nodes.push(new void(frame f, transform t, transform T, pair m, pair M) { + if(srcCopy.nodes.length > 0) { + nodes.push(node(new void(frame f, transform t, transform T, + pair m, pair M) { add(f,srcCopy.fit(t,T*srcCopy.T,m,M),group,filltype,above); - }); + })); + } if(srcCopy.nodes3.length > 0) { - nodes3.push(new void(frame f, transform3 t, transform3 T3, picture pic, - projection P, triple m, triple M) { + nodes3.push(node3(new void(frame f, transform3 t, transform3 T3, + picture pic, projection P, triple m, triple M) + { add(f,srcCopy.fit3(t,T3*srcCopy.T3,pic,P,m,M),group,above); - }); + })); } legend.append(src.legend); @@ -1042,6 +1073,14 @@ void size(picture pic=currentpicture, real x, real y=x, pic.size(x,y,keepAspect); } +void size(picture pic=currentpicture, transform t) +{ + if(pic.empty3()) { + pair z=size(pic.fit(t)); + pic.size(z.x,z.y); + } +} + void size3(picture pic=currentpicture, real x, real y=x, real z=y, bool keepAspect=pic.keepAspect) { @@ -1408,7 +1447,6 @@ void clip(picture pic=currentpicture, path[] g, bool stroke=false, { if(copy) g=copy(g); - //pic.userClip(min(g),max(g)); pic.clip(min(g), max(g), new void(frame f, transform t) { clip(f,t*g,stroke,fillrule,false); diff --git a/Build/source/utils/asymptote/base/plain_scaling.asy b/Build/source/utils/asymptote/base/plain_scaling.asy index fd2dbae07c0..5bed3338e9b 100644 --- a/Build/source/utils/asymptote/base/plain_scaling.asy +++ b/Build/source/utils/asymptote/base/plain_scaling.asy @@ -192,19 +192,30 @@ real max(real M, scaling s, coord[] c) { return M; } -// Calculate the sizing constants for the given array and maximum size. +import simplex; + +/* + Calculate the sizing constants for the given array and maximum size. + Solve the two-variable linear programming problem using the simplex method. + This problem is specialized in that the second variable, "b", does not have + a non-negativity condition, and the first variable, "a", is the quantity + being maximized. +*/ real calculateScaling(string dir, coord[] m, coord[] M, real size, bool warn=true) { - access simplex; - simplex.problem p=new simplex.problem; - + real[][] A; + real[] b; + real[] c=new real[] {-1,0,0}; + void addMinCoord(coord c) { // (a*user + b) + truesize >= 0: - p.addRestriction(c.user,1,c.truesize); - } + A.push(new real[] {c.user,1,-1}); + b.push(-c.truesize); + } void addMaxCoord(coord c) { // (a*user + b) + truesize <= size: - p.addRestriction(-c.user,-1,size-c.truesize); + A.push(new real[] {-c.user,-1,1}); + b.push(c.truesize-size); } for (int i=0; i < m.length; ++i) @@ -212,11 +223,12 @@ real calculateScaling(string dir, coord[] m, coord[] M, real size, for (int i=0; i < M.length; ++i) addMaxCoord(M[i]); - int status=p.optimize(); - if(status == simplex.problem.OPTIMAL) { - // TODO: Could just be return a; - return scaling.build(p.a(),p.b()).a; - } else if(status == simplex.problem.UNBOUNDED) { + int[] s=array(A.length,1); + simplex S=simplex(c,A,s,b); + + if(S.case == S.OPTIMAL) { + return S.x[0]; + } else if(S.case == S.UNBOUNDED) { if(warn) warning("unbounded",dir+" scaling in picture unbounded"); return 0; } else { diff --git a/Build/source/utils/asymptote/base/plain_shipout.asy b/Build/source/utils/asymptote/base/plain_shipout.asy index 972f44963c3..4562205966d 100644 --- a/Build/source/utils/asymptote/base/plain_shipout.asy +++ b/Build/source/utils/asymptote/base/plain_shipout.asy @@ -1,6 +1,10 @@ // Default file prefix used for inline LaTeX mode string defaultfilename; +file _outpipe; +if(settings.xasy) + _outpipe=output(mode="pipe"); + string[] file3; string outprefix(string prefix=defaultfilename) { @@ -14,8 +18,6 @@ string outformat(string format="") return format; } -bool shipped; // Was a picture or frame already shipped out? - frame currentpatterns; frame Portrait(frame f) {return f;}; @@ -58,24 +60,37 @@ frame enclose(string prefix=defaultfilename, object F, string format="") } return F.f; } -include plain_xasy; +void deconstruct(picture pic=currentpicture) +{ + frame f; + transform t=pic.calculateTransform(); + if(currentpicture.fitter == null) + f=pic.fit(t); + else + f=pic.fit(); + deconstruct(f,currentpatterns,t); +} + +bool implicitshipout=false; void shipout(string prefix=defaultfilename, frame f, string format="", bool wait=false, bool view=true, string options="", string script="", - light light=currentlight, projection P=currentprojection) + light light=currentlight, projection P=currentprojection, + transform t=identity) { if(is3D(f)) { f=enclose(prefix,embed3(prefix,f,format,options,script,light,P)); if(settings.render != 0 && !prc(format)) { - shipped=true; return; } } - if(inXasyMode) { - erase(); - add(f,group=false); + if(settings.xasy || (!implicitshipout && prefix == defaultfilename)) { + if(prefix == defaultfilename) { + currentpicture.clear(); + add(f,group=false); + } return; } @@ -84,9 +99,7 @@ void shipout(string prefix=defaultfilename, frame f, int limit=2000; if(abs(m.x) > limit || abs(m.y) > limit) f=shift(-m)*f; - shipout(prefix,f,currentpatterns,format,wait,view, - xformStack.empty() ? null : xformStack.pop0); - shipped=true; + _shipout(prefix,f,currentpatterns,format,wait,view,t); } void shipout(string prefix=defaultfilename, picture pic=currentpicture, @@ -95,6 +108,7 @@ void shipout(string prefix=defaultfilename, picture pic=currentpicture, string options="", string script="", light light=currentlight, projection P=currentprojection) { + pic.uptodate=true; if(!uptodate()) { bool inlinetex=settings.inlinetex; bool prc=prc(format); @@ -110,14 +124,17 @@ void shipout(string prefix=defaultfilename, picture pic=currentpicture, } settings.inlinetex=settings.inlineimage; } - frame f=pic.fit(prefix,format,view=view,options,script,light,P); + frame f; + transform t=pic.calculateTransform(); + if(currentpicture.fitter == null) + f=pic.fit(t); + else + f=pic.fit(prefix,format,view=view,options,script,light,P); + if(!prconly() && (!pic.empty2() || settings.render == 0 || prc || empty3)) - shipout(prefix,orientation(f),format,wait,view); + shipout(prefix,orientation(f),format,wait,view,t); settings.inlinetex=inlinetex; } - - pic.uptodate=true; - shipped=true; } void newpage(picture pic=currentpicture) diff --git a/Build/source/utils/asymptote/base/plain_xasy.asy b/Build/source/utils/asymptote/base/plain_xasy.asy deleted file mode 100644 index efc739d8385..00000000000 --- a/Build/source/utils/asymptote/base/plain_xasy.asy +++ /dev/null @@ -1,132 +0,0 @@ -restricted bool inXasyMode=false; -bool diagnostics=false; -void report(string text) -{ - if(diagnostics) - write(text); -} -void report(transform t) -{ - if(diagnostics) - write(t); -} -void report(int i) -{ - if(diagnostics) - write(i); -} -void initXasyMode() -{ - size(0,0); - inXasyMode=true; -} - -void exitXasyMode() -{ - inXasyMode=false; -} -private picture[] tempStore; -private picture newPic; -void startScript() -{ - tempStore.push(currentpicture.copy()); - newPic=new picture; - currentpicture=newPic; -} - -void endScript() -{ - if(tempStore.length < 1) { - abort("endScript() without matching beginScript()"); - } else { - currentpicture=tempStore.pop(); - add(currentpicture,newPic.fit(),group=false); - } - shipped=false; -} - -struct indexedTransform { - int index; - transform t; - bool active; - void operator init(int index, transform t, bool active=true) { - this.index=index; - this.t=t; - this.active=active; - } -} - -struct framedTransformStack { - struct transact { - transform t; - bool active; - void operator init(transform t, bool active=true) { - this.t=t; - this.active=active; - } - void operator init(indexedTransform i){ - this.t=i.t; - this.active=i.active; - } - void operator init() { - this.t=identity(); - this.active=true; - } - } - private transact[] stack; - private int[] frames; - private int stackBase=0; - transform pop() { - if(stack.length == 0) - return identity(); - else { - transform popped=stack[0].t; - stack.delete(0); - report("Popped"); - report(popped); - return popped; - } - } - - transform pop0() { - if(stack.length == 0) - return identity(); - else { - static transform zerotransform=(0,0,0,0,0,0); - transform popped=stack[0].active ? stack[0].t : zerotransform; - stack.delete(0); - report("Popped"); - report(popped); - return popped; - } - } - - void push(transform t, bool Active=true) { - report("Pushed"); - report(t); - stack.push(transact(t,Active)); - } - - void add(... indexedTransform[] tList) { - transact[] toPush; - for(int a=0; a < tList.length; ++a) - toPush[tList[a].index]=transact(tList[a]); - for(int a=0; a < toPush.length; ++a) - if(!toPush.initialized(a)) - toPush[a]=transact(); - report("Added"); - report(toPush.length); - stack.append(toPush); - } - - bool empty() { - return stack.length == 0; - } -} - -framedTransformStack xformStack; - -void deconstruct(picture pic=currentpicture, real magnification=1) -{ - deconstruct(pic.fit(),currentpatterns,magnification,xformStack.pop); -} diff --git a/Build/source/utils/asymptote/base/rational.asy b/Build/source/utils/asymptote/base/rational.asy index 0fa4ee36afa..a142063806d 100644 --- a/Build/source/utils/asymptote/base/rational.asy +++ b/Build/source/utils/asymptote/base/rational.asy @@ -87,6 +87,11 @@ bool operator ==(rational r, rational s) return r.p == s.p && r.q == s.q; } +bool operator !=(rational r, rational s) +{ + return r.p != s.p || r.q != s.q; +} + bool operator <(rational r, rational s) { return r.p*s.q-s.p*r.q < 0; diff --git a/Build/source/utils/asymptote/base/rationalSimplex.asy b/Build/source/utils/asymptote/base/rationalSimplex.asy new file mode 100644 index 00000000000..fd7128b8025 --- /dev/null +++ b/Build/source/utils/asymptote/base/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); +*/ diff --git a/Build/source/utils/asymptote/base/simplex.asy b/Build/source/utils/asymptote/base/simplex.asy index 43920109ff4..16a495a082e 100644 --- a/Build/source/utils/asymptote/base/simplex.asy +++ b/Build/source/utils/asymptote/base/simplex.asy @@ -1,274 +1,253 @@ -/***** - * simplex.asy - * Andy Hammerlindl 2004/07/27 - * - * Solves the two-variable linear programming problem using the simplex method. - * This problem is specialized in that the second variable, "b", does not have - * a non-negativity condition, and the first variable, "a", is the quantity - * being maximized. - * Correct execution of the algorithm also assumes that the coefficient of "b" - * will be +1 or -1 in every added restriction, and that the problem can be - * initialized to a valid state by pivoting b with one of the slack - * variables. This assumption may in fact be incorrect. - *****/ - -private real infinity=sqrt(0.25*realMax); - -struct problem { - typedef int var; - static var VAR_A = 0; - static var VAR_B = 1; - - static int OPTIMAL = -1; - static var UNBOUNDED = -2; - static int INVALID = -3; - - struct row { - real c, t[]; - } - - // The variables of the rows. - // Initialized for the two variable problem. - var[] v = {VAR_A, VAR_B}; - - // The rows of equalities. - row rowA() { - row r = new row; - r.c = 0; - r.t = new real[] {1, 0}; - return r; - } - row rowB() { - row r = new row; - r.c = 0; - r.t = new real[] {0, 1}; - return r; - } - row[] rows = {rowA(), rowB()}; - - // The number of original variables. - int n = rows.length; - - // Pivot the variable v[col] with vp. - void pivot(int col, var vp) - { - var vc=v[col]; - - // Recalculate rows v[col] and vp for the pivot-swap. - row rvc = rows[vc], rvp = rows[vp]; - real factor=1/rvp.t[col]; // NOTE: Handle rvp.t[col] == 0 case. - rvc.c=-rvp.c*factor; - rvp.c=0; - rvc.t=-rvp.t*factor; - rvp.t *= 0; - rvc.t[col]=factor; - rvp.t[col]=1; - - var a=min(vc,vp); - var b=max(vc,vp); - - // Recalculate the rows other than the two used for the above pivot. - for (var i = 0; i < a; ++i) { - row r=rows[i]; - real m = r.t[col]; - r.c += m*rvc.c; - r.t += m*rvc.t; - r.t[col]=m*factor; +// General simplex solver written by John C. Bowman and Pouria Ramazi, 2018. + +struct simplex { + static int OPTIMAL=0; + static int UNBOUNDED=1; + static int INFEASIBLE=2; + + int case; + real[] x; + real cost; + + int m,n; + int J; + real epsilonA; + + // Row reduce based on pivot E[I][J] + void rowreduce(real[][] E, int N, int I, int J) { + real[] EI=E[I]; + real v=EI[J]; + for(int j=0; j < J; ++j) EI[j] /= v; + EI[J]=1.0; + for(int j=J+1; j <= N; ++j) EI[j] /= v; + + for(int i=0; i < I; ++i) { + real[] Ei=E[i]; + real EiJ=Ei[J]; + for(int j=0; j < J; ++j) + Ei[j] -= EI[j]*EiJ; + Ei[J]=0.0; + for(int j=J+1; j <= N; ++j) + Ei[j] -= EI[j]*EiJ; } - for (var i = a+1; i < b; ++i) { - row r=rows[i]; - real m = r.t[col]; - r.c += m*rvc.c; - r.t += m*rvc.t; - r.t[col]=m*factor; + for(int i=I+1; i <= m; ++i) { + real[] Ei=E[i]; + real EiJ=Ei[J]; + for(int j=0; j < J; ++j) + Ei[j] -= EI[j]*EiJ; + Ei[J]=0.0; + for(int j=J+1; j <= N; ++j) + Ei[j] -= EI[j]*EiJ; } - for (var i = b+1; i < rows.length; ++i) { - row r=rows[i]; - real m = r.t[col]; - r.c += m*rvc.c; - r.t += m*rvc.t; - r.t[col]=m*factor; - } - - // Relabel the vars. - v[col] = vp; } - // As b does not have a non-negativity condition, it must initially be - // pivoted out for a variable that does. This selects the initial - // variable to pivot with b. It also assumes that there is a valid - // solution with a == 0 to the linear programming problem, and if so, it - // picks a pivot to get to that state. In our case, a == 0 corresponds to - // a picture with the user coordinates shrunk down to zero, and if that - // doesn't fit, nothing will. - var initVar() - { - real min=infinity, max=-infinity; - var argmin=0, argmax=0; - - for (var i = 2; i < rows.length; ++i) { - row r=rows[i]; - if (r.t[VAR_B] > 0) { - real val=r.c/r.t[VAR_B]; - if (val < min) { - min=val; - argmin=i; + int iterate(real[][] E, int N, int[] Bindices) { + while(true) { + // Find first negative entry in bottom (reduced cost) row + real[] Em=E[m]; + for(J=0; J < N; ++J) + if(Em[J] < 0) break; + + if(J == N) + return 0; + + 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; + I=i; + break; } - } else if (r.t[VAR_B] < 0) { - real val=r.c/r.t[VAR_B]; - if (val > max) { - max=val; - argmax=i; + } + for(int i=I+1; i < m; ++i) { + real e=E[i][J]; + if(e > epsilonA) { + real v=E[i][N]/e; + if(v <= M) {M=v; I=i;} } } - } + if(I == -1) + return UNBOUNDED; // Can only happen in Phase 2. - // If b has a minimal value, choose a pivot that will give b its minimal - // value. Otherwise, if b has maximal value, choose a pivot to give b its - // maximal value. - return argmin != 0 ? argmin : - argmax != 0 ? argmax : - UNBOUNDED; - } + Bindices[I]=J; - // Initialize the linear program problem by moving into an acceptable state - // this assumes that b is unrestrained and is the second variable. - // NOTE: Works in limited cases, may be bug-ridden. - void init() - { - // Find the lowest constant term in the equations. - var lowest = 0; - for (var i = 2; i < rows.length; ++i) { - if (rows[i].c < rows[lowest].c) - lowest = i; + // Generate new tableau + rowreduce(E,N,I,J); } - - // Pivot if necessary. - if (lowest != 0) - pivot(VAR_B, lowest); + return 0; } - // Selects a column to pivot on. Returns OPTIMAL if the current state is - // optimal. Assumes we are optimizing the first row. - int selectColumn() - { - int i=find(rows[0].t > 0,1); - return (i >= 0) ? i : OPTIMAL; - } - - // Select the new variable associated with a pivot on the column given. - // Returns UNBOUNDED if the space is unbounded. - var selectVar(int col) - { - // We assume that the first two vars (a and b) once swapped out, won't be - // swapped back in. This finds the variable which gives the tightest - // non-negativity condition restricting our optimization. This turns - // out to be the max of c/t[col]. Note that as c is positive, and - // t[col] is negative, all c/t[col] will be negative, so we are finding - // the smallest in magnitude. - var vp=UNBOUNDED; - real max=-infinity; - for (int i = 2; i < rows.length; ++i) { - row r=rows[i]; - if(r.c < max*r.t[col]) { - max=r.c/r.t[col]; vp=i; + // 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) { + static real epsilon=sqrt(realEpsilon); + epsilonA=epsilon*norm(A); + + // 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; + real[][] E=new real[m+1][N+1]; + real[] Em=E[m]; + + for(int j=0; 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) { + for(int j=0; j < n; ++j) { + real Aij=Ai[j]; + Ei[j]=Aij; + Em[j] -= Aij; + } + } else { + for(int j=0; j < n; ++j) { + real Aij=-Ai[j]; + Ei[j]=Aij; + Em[j] -= Aij; + } } } - - return vp; - } - // Checks that the rows are in a valid state. - bool valid() - { - // Checks that constants are valid. - bool validConstants() { - for (int i = 0; i < rows.length; ++i) - // Do not test the row for b, as it does not have a non-negativity - // condition. - if (i != VAR_B && rows[i].c < 0) - return false; - return true; + 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; + } } - // Check a variable to see if its row is simple. - // NOTE: Simple rows could be optimized out, since they are not really - // used. - bool validVar(int col) { - - var vc = v[col]; - row rvc = rows[vc]; - - if (rvc.c != 0) - return false; - for (int i = 0; i < n; ++i) - if (rvc.t[i] != (i == col ? 1 : 0)) - return false; - - return true; + real sum=0; + for(int i=0; i < m; ++i) { + real B=abs(b[i]); + E[i][N]=B; + sum -= B; } - - if (!validConstants()) { - return false; + Em[N]=sum; + + if(phase1) + for(int j=0; j < m; ++j) + Em[n+j]=0.0; + + int[] Bindices=sequence(new int(int x){return x;},m)+n; + + if(phase1) { + iterate(E,N,Bindices); + + if(abs(Em[J]) > epsilonA) { + case=INFEASIBLE; + return; + } } - for (int i = 0; i < n; ++i) - if (!validVar(i)) { - return false; + + 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) { + 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]; + real[] Dip=D[ip]; + real[] Ei=E[i]; + for(int j=0; j < n; ++j) + Dip[j]=Ei[j]; + Dip[n]=Ei[N]; + ++ip; } - return true; - } - + real[] Dip=D[ip]; + real[] Em=E[m]; + for(int j=0; j < n; ++j) + Dip[j]=Em[j]; + Dip[n]=Em[N]; - // Perform the algorithm to find the optimal solution. Returns OPTIMAL, - // UNBOUNDED, or INVALID (if no solution is possible). - int optimize() - { - // Put into a valid state to begin and pivot b out. - var iv=initVar(); - if (iv == UNBOUNDED) - return iv; - pivot(VAR_B, iv); + m=ip; - if (!valid()) - return INVALID; + 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; + } - while(true) { - int col = selectColumn(); - - if (col == OPTIMAL) - return col; - var vp = selectVar(col); - - if (vp == UNBOUNDED) - return vp; - - pivot(col, vp); + // Done with Phase 1 + } + + real 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; } - // Shouldn't reach here. - return INVALID; - } + for(int j=0; j < n; ++j) + x[j]=0; - // Add a restriction to the problem: - // t1*a + t2*b + c >= 0 - void addRestriction(real t1, real t2, real c) - { - row r = new row; - r.c = c; - r.t = new real[] {t1, t2}; - rows.push(r); - } + for(int k=0; k < m; ++k) + x[Bindices[k]]=D[k][n]; - // Return the value of a computed. - real a() - { - return rows[VAR_A].c; + cost=-Dm[n]; + case=OPTIMAL; } - // Return the value of b computed. - real b() - { - return rows[VAR_B].c; + // 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 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, int[] s, real[] 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) + if(s[i] != 0) ++count; + + real[][] a=new real[m][n+count]; + + for(int i=0; i < m; ++i) { + real[] ai=a[i]; + real[] Ai=A[i]; + for(int j=0; j < n; ++j) { + ai[j]=Ai[j]; + } + } + + int k=0; + + for(int i=0; i < m; ++i) { + real[] 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); // TODO: Check + bool phase1=true; + operator init(concat(c,array(count,0.0)),a,b,phase1); + + if(case == OPTIMAL) + x.delete(n,n+count-1); } } diff --git a/Build/source/utils/asymptote/base/slide.asy b/Build/source/utils/asymptote/base/slide.asy index 6195a97ec01..8a787c37a5b 100644 --- a/Build/source/utils/asymptote/base/slide.asy +++ b/Build/source/utils/asymptote/base/slide.asy @@ -225,7 +225,7 @@ void erasestep(int erasenode) { for(int i=0; i < firstnode.length; ++i) { for(int j=firstnode[i]; j <= lastnode[i]; ++j) { tex(bulletcolor(oldbulletcolor)); - currentpicture.add(currentpicture.nodes[j]); + currentpicture.add(currentpicture.nodes[j].d); } } firstnode.push(currentpicture.nodes.length-1); diff --git a/Build/source/utils/asymptote/base/three.asy b/Build/source/utils/asymptote/base/three.asy index 01eec97976f..742ceda1c04 100644 --- a/Build/source/utils/asymptote/base/three.asy +++ b/Build/source/utils/asymptote/base/three.asy @@ -1,6 +1,7 @@ private import math; -if(inXasyMode) settings.render=0; +if(settings.xasy) + settings.render=0; if(prc0()) { if(!latex()) settings.prc=false; @@ -2916,9 +2917,9 @@ object embed(string prefix=outprefix(), string label=prefix, if((preview || (prc && settings.render == 0)) && settings.embed) { image=prefix; if(settings.inlinetex) image += "_0"; - if(!preview && !shipped && !S.pic2.empty2()) { + if(!preview && !S.pic2.empty2()) { transform T=S.pic2.scaling(S.width,S.height); - shipout(image,S.pic2.fit(T),newframe,nativeformat(),false,false,null); + _shipout(image,S.pic2.fit(T),newframe,nativeformat(),false,false); } image += "."+nativeformat(); -- cgit v1.2.3