// A linear scale, with optional autoscaling of minimum and maximum values,
// scaling factor s and intercept.
scaleT Linear(bool automin=false, bool automax=automin, real s=1,
real intercept=0)
{
real sinv=1/s;
scalefcn T,Tinv;
if(s == 1 && intercept == 0)
T=Tinv=identity;
else {
T=new real(real x) {return (x-intercept)*s;};
Tinv=new real(real x) {return x*sinv+intercept;};
}
return scaleT(T,Tinv,logarithmic=false,automin,automax);
}
// A logarithmic scale, with optional autoscaling of minimum and maximum
// values.
scaleT Log(bool automin=false, bool automax=automin)
{
return scaleT(Log.T,Log.Tinv,logarithmic=true,automin,automax);
}
// A "broken" linear axis omitting the segment [a,b].
scaleT Broken(real a, real b, bool automin=false, bool automax=automin)
{
real skip=b-a;
real T(real x) {
if(x <= a) return x;
if(x <= b) return a;
return x-skip;
}
real Tinv(real x) {
if(x <= a) return x;
return x+skip;
}
return scaleT(T,Tinv,logarithmic=false,automin,automax);
}
// A "broken" logarithmic axis omitting the segment [a,b], where a and b are
// automatically rounded to the nearest integral power of the base.
scaleT BrokenLog(real a, real b, bool automin=false, bool automax=automin)
{
real A=round(Log.T(a));
real B=round(Log.T(b));
a=Log.Tinv(A);
b=Log.Tinv(B);
real skip=B-A;
real T(real x) {
if(x <= a) return Log.T(x);
if(x <= b) return A;
return Log.T(x)-skip;
}
real Tinv(real x) {
real X=Log.Tinv(x);
if(X <= a) return X;
return Log.Tinv(x+skip);
}
return scaleT(T,Tinv,logarithmic=true,automin,automax);
}
struct scientific
{
int sign;
real mantissa;
int exponent;
int ceil() {return sign*ceil(mantissa);}
real scale(real x, real exp) {
static real max=0.1*realMax;
static real limit=-log10(max);
return x*(exp > limit ? 10^-exp : max);
}
real ceil(real x, real exp) {return ceil(sign*scale(abs(x),exp));}
real floor(real x, real exp) {return floor(sign*scale(abs(x),exp));}
}
// Convert x to scientific notation
scientific scientific(real x)
{
scientific s;
s.sign=sgn(x);
x=abs(x);
if(x == 0) {s.mantissa=0; s.exponent=-intMax; return s;}
real logx=log10(x);
s.exponent=floor(logx);
s.mantissa=s.scale(x,s.exponent);
return s;
}
// Autoscale limits and tick divisor.
struct bounds {
real min;
real max;
// Compute tick divisors.
int[] divisors(int a, int b)
{
int[] dlist;
int n=b-a;
dlist[0]=1;
if(n == 1) {dlist[1]=10; dlist[2]=100; return dlist;}
if(n == 2) {dlist[1]=2; return dlist;}
int sqrtn=floor(sqrt(n));
int i=0;
for(int d=2; d <= sqrtn; ++d)
if(n % d == 0 && (a*b >= 0 || b % (n/d) == 0)) dlist[++i]=d;
for(int d=sqrtn; d >= 1; --d)
if(n % d == 0 && (a*b >= 0 || b % d == 0)) dlist[++i]=quotient(n,d);
return dlist;
}
real upscale(real b, real a)
{
if(b <= 5) b=5;
else if (b > 10 && a >= 0 && b <= 12) b=12;
else if (b > 10 && (a >= 0 || 15 % -a == 0) && b <= 15) b=15;
else b=ceil(b/10)*10;
return b;
}
// Format tick values as integral powers of base; otherwise with DefaultFormat.
ticklabel DefaultLogFormat(int base) {
return new string(real x) {
string exponent=format("%.4f",log(x)/log(base));
return find(exponent,".") == -1 ?
"$"+(string) base+"^{"+exponent+"}$" : format(x);
};
}
// Format all tick values as powers of base.
ticklabel LogFormat(int base) {
return new string(real x) {
return format("$"+(string) base+"^{%g}$",log(x)/log(base));
};
}
// The default direction specifier.
pair zero(real) {return 0;}
struct ticklocate {
real a,b; // Tick values at point(g,0), point(g,length(g)).
autoscaleT S; // Autoscaling transformation.
pair dir(real t); // Absolute 2D tick direction.
triple dir3(real t); // Absolute 3D tick direction.
real time(real v); // Returns the time corresponding to the value v.
ticklocate copy() {
ticklocate T=new ticklocate;
T.a=a;
T.b=b;
T.S=S.copy();
T.dir=dir;
T.dir3=dir3;
T.time=time;
return T;
}
}
autoscaleT defaultS;
typedef real valuetime(real);
valuetime linear(scalefcn S=identity, real Min, real Max)
{
real factor=Max == Min ? 0.0 : 1.0/(Max-Min);
return new real(real v) {return (S(v)-Min)*factor;};
}
ticklocate ticklocate(real a, real b, autoscaleT S=defaultS,
real tickmin=-infinity, real tickmax=infinity,
real time(real)=null, pair dir(real)=zero)
{
if((valuetime) time == null) time=linear(S.T(),a,b);
ticklocate locate;
locate.a=a;
locate.b=b;
locate.S=S.copy();
if(finite(tickmin)) locate.S.tickMin=tickmin;
if(finite(tickmax)) locate.S.tickMax=tickmax;
locate.time=time;
locate.dir=dir;
return locate;
}
private struct locateT {
real t; // tick location time
pair Z; // tick location in frame coordinates
pair pathdir; // path direction in frame coordinates
pair dir; // tick direction in frame coordinates
void dir(transform T, path g, ticklocate locate, real t) {
pathdir=unit(shiftless(T)*dir(g,t));
pair Dir=locate.dir(t);
dir=Dir == 0 ? -I*pathdir : unit(Dir);
}
// Locate the desired position of a tick along a path.
void calc(transform T, path g, ticklocate locate, real val) {
t=locate.time(val);
Z=T*point(g,t);
dir(T,g,locate,t);
}
}
// Check the tick coverage of a linear axis.
bool axiscoverage(int N, transform T, path g, ticklocate locate, real Step,
pair side, int sign, real Size, Label F, ticklabel ticklabel,
real norm, real limit)
{
real coverage=0;
bool loop=cyclic(g);
real a=locate.S.Tinv(locate.a);
real b=locate.S.Tinv(locate.b);
real tickmin=finite(locate.S.tickMin) ? locate.S.Tinv(locate.S.tickMin) : a;
if(Size > 0) {
int count=0;
if(loop) count=N+1;
else {
for(int i=0; i <= N; ++i) {
real val=tickmin+i*Step;
if(val >= a && val <= b)
++count;
}
}
if(count > 0) limit /= count;
for(int i=0; i <= N; ++i) {
real val=tickmin+i*Step;
if(loop || (val >= a && val <= b)) {
frame d;
pair dir=labeltick(d,T,g,locate,val,side,sign,Size,ticklabel,F,norm);
if(abs(dot(size(d),dir)) > limit) return false;
}
}
}
return true;
}
// Check the tick coverage of a logarithmic axis.
bool logaxiscoverage(int N, transform T, path g, ticklocate locate, pair side,
int sign, real Size, Label F, ticklabel ticklabel,
real limit, int first, int last)
{
bool loop=cyclic(g);
real coverage=0;
real a=locate.a;
real b=locate.b;
int count=0;
for(int i=first-1; i <= last+1; i += N) {
if(loop || i >= a && i <= b)
++count;
}
if(count > 0) limit /= count;
for(int i=first-1; i <= last+1; i += N) {
if(loop || i >= a && i <= b) {
frame d;
pair dir=labeltick(d,T,g,locate,i,side,sign,Size,ticklabel,F);
if(abs(dot(size(d),dir)) > limit) return false;
}
}
return true;
}
struct tickvalues {
real[] major;
real[] minor;
int N; // For logarithmic axes: number of decades between tick labels.
}
// Determine a format that distinguishes adjacent pairs of ticks, optionally
// adding trailing zeros.
string autoformat(string format="", real norm ... real[] a)
{
bool trailingzero=(format == trailingzero);
bool signedtrailingzero=(format == signedtrailingzero);
if(!trailingzero && !signedtrailingzero && format != "") return format;
for(int i=0; i < A.length-1; ++i) {
real a=A[i];
real b=A[i+1];
// Check if an extra digit of precision should be added.
string fixedformat="%#."+string(n+1)+"f";
string A=format(fixedformat,a);
string B=format(fixedformat,b);
if(substr(A,length(A)-1,1) != "0" || substr(B,length(B)-1,1) != "0") {
a *= 0.1;
b *= 0.1;
}
if(a != b) {
while(format(format,a) == format(format,b))
format=defaultformat(++n,trailing,Fixed);
}
}
if(n == 0) return defaultformat;
return format;
}
// Automatic tick generation routine.
tickvalues generateticks(int sign, Label F="", ticklabel ticklabel=null,
int N, int n=0, real Step=0, real step=0,
real Size=0, real size=0,
transform T, pair side, path g, real limit,
pen p, ticklocate locate, int[] divisor,
bool opposite)
{
tickvalues tickvalues;
sign=opposite ? -sign : sign;
if(Size == 0) Size=Ticksize;
if(size == 0) size=ticksize;
F=F.copy();
F.p(p);
if(Size > 0) {
for(int i=0; i <= N; ++i) {
real val=tickmin+i*Step;
if(val >= a && val <= b)
tickvalues.major.push(val);
if(size > 0 && step > 0) {
real iStep=i*Step;
real jstop=(len-iStep)/step;
for(int j=1; j < n && j <= jstop; ++j) {
real val=tickmin+iStep+j*step;
if(val >= a && val <= b)
tickvalues.minor.push(val);
}
}
}
}
begingroup(f);
if(opposite) draw(f,G,p);
else draw(f,margin(G,p).g,p,arrow);
for(int i=(begin ? 0 : 1); i < (end ? Ticks.length : Ticks.length-1); ++i) {
real val=T(Ticks[i]);
if(val >= a && val <= b)
drawtick(f,t,g,g2,locate,val,Size,sign,pTick,extend);
}
for(int i=0; i < ticks.length; ++i) {
real val=T(ticks[i]);
if(val >= a && val <= b)
drawtick(f,t,g,g2,locate,val,size,sign,ptick,extend);
}
endgroup(f);
if(N == 0) N=1;
if(Size > 0 && !opposite) {
for(int i=(beginlabel ? 0 : 1);
i < (endlabel ? Ticks.length : Ticks.length-1); i += N) {
real val=T(Ticks[i]);
if(val >= a && val <= b) {
ticklabels=true;
labeltick(f,t,g,locate,val,side,sign,Size,ticklabel,F,norm);
}
}
}
if(L.s != "" && !opposite)
labelaxis(f,t,L,G,locate,sign,ticklabels);
};
}
// Optional routine to allow modification of auto-generated tick values.
typedef tickvalues tickmodifier(tickvalues);
tickvalues None(tickvalues v) {return v;}
// Tickmodifier that removes all ticks in the intervals [a[i],b[i]].
tickmodifier OmitTickIntervals(real[] a, real[] b) {
return new tickvalues(tickvalues v) {
if(a.length != b.length) abort(differentlengths);
void omit(real[] A) {
if(A.length != 0) {
real norm=max(abs(A));
for(int i=0; i < a.length; ++i) {
int j;
while((j=find(A > a[i]-zerotickfuzz*norm
& A < b[i]+zerotickfuzz*norm)) >= 0) {
A.delete(j);
}
}
}
}
omit(v.major);
omit(v.minor);
return v;
};
}
// Tickmodifier that removes all ticks in the interval [a,b].
tickmodifier OmitTickInterval(real a, real b) {
return OmitTickIntervals(new real[] {a}, new real[] {b});
}
// Tickmodifier that removes the specified ticks.
tickmodifier OmitTick(... real[] x) {
return OmitTickIntervals(x,x);
}
tickmodifier NoZero=OmitTick(0);
tickmodifier Break(real, real)=OmitTickInterval;
// Automatic tick construction routine.
ticks Ticks(int sign, Label F="", ticklabel ticklabel=null,
bool beginlabel=true, bool endlabel=true,
int N, int n=0, real Step=0, real step=0,
bool begin=true, bool end=true, tickmodifier modify=None,
real Size=0, real size=0, bool extend=false,
pen pTick=nullpen, pen ptick=nullpen)
{
return new void(frame f, transform T, Label L, pair side, path g, path g2,
pen p, arrowbar arrow, margin margin, ticklocate locate,
int[] divisor, bool opposite) {
real limit=Step == 0 ? axiscoverage*arclength(T*g) : 0;
tickvalues values=modify(generateticks(sign,F,ticklabel,N,n,Step,step,
Size,size,T,side,g,
limit,p,locate,divisor,opposite));
// Structure used to communicate axis and autoscale settings to tick routines.
struct axisT {
int type; // -1 = min, 0 = given value, 1 = max, 2 = min/max
int type2; // for 3D axis
real value;
real value2;
pair side; // 2D tick label direction relative to path (left or right)
real position; // label position along axis
align align; // default axis label alignment and 3D tick label direction
int[] xdivisor;
int[] ydivisor;
int[] zdivisor;
bool extend; // extend axis to graph boundary?
};
real xtrans(transform t, real x)
{
return (t*(x,0)).x;
}
real ytrans(transform t, real y)
{
return (t*(0,y)).y;
}
// An internal routine to draw an x axis at a particular y value.
void xaxisAt(picture pic=currentpicture, Label L="", axis axis,
real xmin=-infinity, real xmax=infinity, pen p=currentpen,
ticks ticks=NoTicks, arrowbar arrow=None, margin margin=NoMargin,
bool above=true, bool opposite=false)
{
real y=axis.value;
real y2;
Label L=L.copy();
int[] divisor=copy(axis.xdivisor);
pair side=axis.side;
int type=axis.type;
// Process any queued y axis bound calculation requests.
for(int i=0; i < pic.scale.y.bound.length; ++i)
pic.scale.y.bound[i]();
pic.scale.y.bound.delete();
bounds();
// Request another x bounds calculation before final picture scaling.
pic.scale.x.bound.push(bounds);
}
// An internal routine to draw a y axis at a particular x value.
void yaxisAt(picture pic=currentpicture, Label L="", axis axis,
real ymin=-infinity, real ymax=infinity, pen p=currentpen,
ticks ticks=NoTicks, arrowbar arrow=None, margin margin=NoMargin,
bool above=true, bool opposite=false)
{
real x=axis.value;
real x2;
Label L=L.copy();
int[] divisor=copy(axis.ydivisor);
pair side=axis.side;
int type=axis.type;
// Process any queued x axis bound calculation requests.
for(int i=0; i < pic.scale.x.bound.length; ++i)
pic.scale.x.bound[i]();
pic.scale.x.bound.delete();
bounds();
// Request another y bounds calculation before final picture scaling.
pic.scale.y.bound.push(bounds);
}
// Set the x limits of a picture.
void xlimits(picture pic=currentpicture, real min=-infinity, real max=infinity,
bool crop=NoCrop)
{
if(min > max) return;
// Set the y limits of a picture.
void ylimits(picture pic=currentpicture, real min=-infinity, real max=infinity,
bool crop=NoCrop)
{
if(min > max) return;
// Crop a picture to the current user-space picture limits.
void crop(picture pic=currentpicture)
{
xlimits(pic,false);
ylimits(pic,false);
if(pic.userSetx() && pic.userSety())
clip(pic,box(pic.userMin(),pic.userMax()));
}
// Restrict the x and y limits to box(min,max).
void limits(picture pic=currentpicture, pair min, pair max, bool crop=NoCrop)
{
xlimits(pic,min.x,max.x);
ylimits(pic,min.y,max.y);
if(crop && pic.userSetx() && pic.userSety())
clip(pic,box(pic.userMin(),pic.userMax()));
}
// Internal routine to autoscale the user limits of a picture.
void autoscale(picture pic=currentpicture, axis axis)
{
if(!pic.scale.set) {
bounds mx,my;
pic.scale.set=true;
// Draw a yaxis at x.
void xequals(picture pic=currentpicture, Label L="", real x,
bool extend=false, real ymin=-infinity, real ymax=infinity,
pen p=currentpen, ticks ticks=NoTicks,
arrowbar arrow=None, margin margin=NoMargin, bool above=true)
{
yaxis(pic,L,XEquals(x,extend),ymin,ymax,p,ticks,arrow,margin,above);
}
// Draw an xaxis at y.
void yequals(picture pic=currentpicture, Label L="", real y,
bool extend=false, real xmin=-infinity, real xmax=infinity,
pen p=currentpen, ticks ticks=NoTicks,
arrowbar arrow=None, margin margin=NoMargin, bool above=true)
{
xaxis(pic,L,YEquals(y,extend),xmin,xmax,p,ticks,arrow,margin,above);
}
real ScaleX(picture pic=currentpicture, real x)
{
return pic.scale.x.T(x);
}
real ScaleY(picture pic=currentpicture, real y)
{
return pic.scale.y.T(y);
}
// Draw a tick of length size at pair z in direction dir using pen p.
void tick(picture pic=currentpicture, pair z, pair dir, real size=Ticksize,
pen p=currentpen)
{
pair z=Scale(pic,z);
pic.add(new void (frame f, transform t) {
pair tz=t*z;
draw(f,tz--tz+unit(dir)*size,p);
});
pic.addPoint(z,p);
pic.addPoint(z,unit(dir)*size,p);
}
// Put a label on the x axis.
void labelx(picture pic=currentpicture, Label L="", explicit pair z,
align align=S, string format="", pen p=currentpen)
{
label(pic,L,Scale(pic,z),z.x,align,format,p);
}
// Put a label on the y axis.
void labely(picture pic=currentpicture, Label L="", explicit pair z,
align align=W, string format="", pen p=currentpen)
{
label(pic,L,Scale(pic,z),z.y,align,format,p);
}
void labely(picture pic=currentpicture, Label L="", real y,
align align=W, string format="", pen p=currentpen)
{
labely(pic,L,(pic.scale.x.scale.logarithmic ? 1 : 0,y),align,format,p);
}
graph graph(interpolate join)
{
return new guide(pair f(real), real a, real b, int n) {
real width=b-a;
return n == 0 ? join(f(a)) :
join(...sequence(new guide(int i) {return f(a+(i/n)*width);},n+1));
};
}
multigraph graph(interpolate join, bool3 cond(real))
{
return new guide[](pair f(real), real a, real b, int n) {
real width=b-a;
if(n == 0) return new guide[] {join(cond(a) ? f(a) : nullpath)};
guide[] G;
guide[] g;
for(int i=0; i < n+1; ++i) {
real t=a+(i/n)*width;
bool3 b=cond(t);
if(b)
g.push(f(t));
else {
if(g.length > 0) {
G.push(join(...g));
g=new guide[] {};
}
if(b == default)
g.push(f(t));
}
}
if(g.length > 0)
G.push(join(...g));
return G;
};
}
// Connect points in z into segments corresponding to consecutive true elements
// of b using interpolation operator join.
path[] segment(pair[] z, bool[] cond, interpolate join=operator --)
{
checkconditionlength(cond.length,z.length);
int[][] segment=segment(cond);
return sequence(new path(int i) {return join(...z[segment[i]]);},
segment.length);
}
pair polar(real r, real theta)
{
return r*expi(theta);
}
guide polargraph(picture pic=currentpicture, real r(real), real a, real b,
int n=ngraph, interpolate join=operator --)
{
return graph(join)(new pair(real theta) {
return Scale(pic,polar(r(theta),theta));
},a,b,n);
}
// Graph a parametric function with the control points specified to match the
// derivative.
guide graphwithderiv(pair f(real), pair fprime(real), real a, real b,
int n=ngraph # 10) {
guide toreturn;
real segmentlen = (b - a) / n;
pair[] points = new pair[n+1];
pair[] derivs = new pair[n+1];
real derivfactor = segmentlen / 3;
for (int i = 0; i <= n; ++i) {
real t = interp(a, b, i/n);
points[i] = f(t);
derivs[i] = fprime(t) * derivfactor;
}
toreturn = points[0];
for (int i = 1; i <= n; ++i) {
toreturn = toreturn .. controls (points[i-1] + derivs[i-1])
and (points[i] - derivs[i]) .. points[i];
}
return toreturn;
}
void errorbar(picture pic, pair z, pair dp, pair dm, pen p=currentpen,
real size=0)
{
real dmx=-abs(dm.x);
real dmy=-abs(dm.y);
real dpx=abs(dp.x);
real dpy=abs(dp.y);
if(dmx != dpx) draw(pic,Scale(pic,z+(dmx,0))--Scale(pic,z+(dpx,0)),p,
Bars(size));
if(dmy != dpy) draw(pic,Scale(pic,z+(0,dmy))--Scale(pic,z+(0,dpy)),p,
Bars(size));
}
void errorbars(picture pic=currentpicture, real[] x, real[] y,
real[] dpy, bool[] cond={}, pen p=currentpen, real size=0)
{
errorbars(pic,x,y,0*x,dpy,cond,p,size);
}
// Return a vector field on path g, specifying the vector as a function of the
// relative position along path g in [0,1].
picture vectorfield(path vector(real), path g, int n, bool truesize=false,
pen p=currentpen, arrowbar arrow=Arrow,
margin margin=PenMargin)
{
picture pic;
for(int i=0; i < n; ++i) {
real x=(n == 1) ? 0.5 : i/(n-1);
if(truesize)
draw(relpoint(g,x),pic,vector(x),p,arrow);
else
draw(pic,shift(relpoint(g,x))*vector(x),p,arrow,margin);
}
return pic;
}
// return a vector field over box(a,b).
picture vectorfield(path vector(pair), pair a, pair b,
int nx=nmesh, int ny=nx, bool truesize=false,
bool cond(pair z)=null, pen p=currentpen,
arrowbar arrow=Arrow, margin margin=PenMargin)
{
picture pic;
real dx=(b.x-a.x)/(nx-1);
real dy=(b.y-a.y)/(ny-1);
bool all=cond == null;
real scale;
