summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/base
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/utils/asymptote/base')
-rw-r--r--Build/source/utils/asymptote/base/graph.asy5
-rw-r--r--Build/source/utils/asymptote/base/graph3.asy131
-rw-r--r--Build/source/utils/asymptote/base/lmfit.asy838
-rw-r--r--Build/source/utils/asymptote/base/math.asy29
-rw-r--r--Build/source/utils/asymptote/base/ode.asy257
-rw-r--r--Build/source/utils/asymptote/base/plain.asy2
-rw-r--r--Build/source/utils/asymptote/base/plain_strings.asy4
-rw-r--r--Build/source/utils/asymptote/base/slide.asy1
-rw-r--r--Build/source/utils/asymptote/base/splinetype.asy4
-rw-r--r--Build/source/utils/asymptote/base/three.asy47
-rw-r--r--Build/source/utils/asymptote/base/three_surface.asy6
11 files changed, 1231 insertions, 93 deletions
diff --git a/Build/source/utils/asymptote/base/graph.asy b/Build/source/utils/asymptote/base/graph.asy
index 805cf15f612..2983137b612 100644
--- a/Build/source/utils/asymptote/base/graph.asy
+++ b/Build/source/utils/asymptote/base/graph.asy
@@ -518,16 +518,17 @@ string autoformat(string format="", real norm ... real[] a)
bool signchange=(A.length > 1 && A[0] < 0 && A[A.length-1] >= 0);
- for(int i=0; i < A.length-1; ++i)
+ for(int i=0; i < A.length; ++i)
if(a[i] < zerotickfuzz*norm) A[i]=a[i]=0;
int n=0;
bool Fixed=find(a >= 1e4-epsilon | (a > 0 & a <= 1e-4-epsilon)) < 0;
+
string Format=defaultformat(4,fixed=Fixed);
if(Fixed && n < 4) {
- for(int i=0; i < A.length-1; ++i) {
+ for(int i=0; i < A.length; ++i) {
real a=A[i];
while(format(defaultformat(n,fixed=Fixed),a) != format(Format,a))
++n;
diff --git a/Build/source/utils/asymptote/base/graph3.asy b/Build/source/utils/asymptote/base/graph3.asy
index a1c8258638f..8fe6c9ed5c6 100644
--- a/Build/source/utils/asymptote/base/graph3.asy
+++ b/Build/source/utils/asymptote/base/graph3.asy
@@ -1591,15 +1591,11 @@ private surface bispline(real[][] z, real[][] p, real[][] q, real[][] r,
real yj=y[j];
real yp=y[j+1];
if(all || condi[j]) {
- triple[][] P={
- {O,O,O,O},
- {O,O,O,O},
- {O,O,O,O},
- {O,O,O,O}};
+ triple[][] P=array(4,array(4,O));
real hy=(yp-yj)/3;
real hxy=hx*hy;
// first x and y directions
- for(int k=0 ; k < 4 ; ++k) {
+ for(int k=0; k < 4; ++k) {
P[k][0] += xi*X;
P[0][k] += yj*Y;
P[k][1] += (xp+2*xi)/3*X;
@@ -1636,34 +1632,38 @@ private surface bispline(real[][] z, real[][] p, real[][] q, real[][] r,
}
// return the surface described by a real matrix f, interpolated with
-// splinetype.
+// xsplinetype and ysplinetype.
surface surface(real[][] f, real[] x, real[] y,
- splinetype splinetype=null, bool[][] cond={})
+ splinetype xsplinetype=null, splinetype ysplinetype=xsplinetype,
+ bool[][] cond={})
{
- if(splinetype == null)
- splinetype=(x[0] == x[x.length-1] && y[0] == y[y.length-1]) ?
- periodic : notaknot;
+ real epsilon=sqrtEpsilon*max(abs(y));
+ if(xsplinetype == null)
+ xsplinetype=(abs(x[0]-x[x.length-1]) <= epsilon) ? periodic : notaknot;
+ if(ysplinetype == null)
+ ysplinetype=(abs(y[0]-y[y.length-1]) <= epsilon) ? periodic : notaknot;
int n=x.length; int m=y.length;
real[][] ft=transpose(f);
real[][] tp=new real[m][];
- for(int j=0; j < m ; ++j)
- tp[j]=splinetype(x,ft[j]);
+ for(int j=0; j < m; ++j)
+ tp[j]=xsplinetype(x,ft[j]);
real[][] q=new real[n][];
- for(int i=0; i < n ; ++i)
- q[i]=splinetype(y,f[i]);
+ for(int i=0; i < n; ++i)
+ q[i]=ysplinetype(y,f[i]);
real[][] qt=transpose(q);
- real[] d1=splinetype(x,qt[0]);
- real[] d2=splinetype(x,qt[m-1]);
+ real[] d1=xsplinetype(x,qt[0]);
+ real[] d2=xsplinetype(x,qt[m-1]);
real[][] r=new real[n][];
- for(int i=0; i < n ; ++i)
- r[i]=clamped(d1[i],d2[i])(y,f[i]);
- return bispline(f,transpose(tp),q,r,x,y,cond);
+ real[][] p=transpose(tp);
+ for(int i=0; i < n; ++i)
+ r[i]=clamped(d1[i],d2[i])(y,p[i]);
+ return bispline(f,p,q,r,x,y,cond);
}
// return the surface described by a real matrix f, interpolated with
-// splinetype.
-surface surface(real[][] f, pair a, pair b, splinetype splinetype,
- bool[][] cond={})
+// xsplinetype and ysplinetype.
+surface surface(real[][] f, pair a, pair b, splinetype xsplinetype,
+ splinetype ysplinetype=xsplinetype, bool[][] cond={})
{
if(!rectangular(f)) abort("matrix is not rectangular");
@@ -1674,7 +1674,7 @@ surface surface(real[][] f, pair a, pair b, splinetype splinetype,
real[] x=uniform(a.x,b.x,nx);
real[] y=uniform(a.y,b.y,ny);
- return surface(f,x,y,splinetype,cond);
+ return surface(f,x,y,xsplinetype,ysplinetype,cond);
}
// return the surface described by a real matrix f, interpolated linearly.
@@ -1732,6 +1732,82 @@ surface surface(triple f(pair z), pair a, pair b, int nu=nmesh, int nv=nu,
return surface(v,active);
}
+// return the surface described by a parametric function f over box(a,b),
+// interpolated with usplinetype and vsplinetype.
+surface surface(triple f(pair z), pair a, pair b, int nu=nmesh, int nv=nu,
+ splinetype[] usplinetype, splinetype[] vsplinetype=Spline,
+ bool cond(pair z)=null)
+{
+ real[] upt=uniform(a.x,b.x,nu);
+ real[] vpt=uniform(a.y,b.y,nv);
+ real[] ipt=sequence(nu+1);
+ real[] jpt=sequence(nv+1);
+ real[][] fx=new real[nu+1][nv+1];
+ real[][] fy=new real[nu+1][nv+1];
+ real[][] fz=new real[nu+1][nv+1];
+
+ bool[][] active;
+ bool all=cond == null;
+ if(!all) active=new bool[nu+1][nv+1];
+
+ real norm;
+ for(int i=0; i <= nu; ++i) {
+ real upti=upt[i];
+ real[] fxi=fx[i];
+ real[] fyi=fy[i];
+ real[] fzi=fz[i];
+ bool[] activei=all ? null : active[i];
+ for(int j=0; j <= nv; ++j) {
+ pair z=(upti,vpt[j]);
+ triple f=(all || (activei[j]=cond(z))) ? f(z) : O;
+ norm=max(norm,abs(f));
+ fxi[j]=f.x;
+ fyi[j]=f.y;
+ fzi[j]=f.z;
+ }
+ }
+
+ real epsilon=sqrtEpsilon*norm;
+
+ if(usplinetype.length == 0) {
+ bool uperiodic(real[][] a) {
+ for(int j=0; j < nv; ++j)
+ if(abs(a[0][j]-a[nu][j]) > epsilon) return false;
+ return true;
+ }
+ usplinetype=new splinetype[] {uperiodic(fx) ? periodic : notaknot,
+ uperiodic(fy) ? periodic : notaknot,
+ uperiodic(fz) ? periodic : notaknot};
+ } else if(usplinetype.length != 3) abort("usplinetype must have length 3");
+
+ if(vsplinetype.length == 0) {
+ bool vperiodic(real[][] a) {
+ for(int i=0; i < nu; ++i)
+ if(abs(a[i][0]-a[i][nv]) > epsilon) return false;
+ return true;
+ }
+ vsplinetype=new splinetype[] {vperiodic(fx) ? periodic : notaknot,
+ vperiodic(fy) ? periodic : notaknot,
+ vperiodic(fz) ? periodic : notaknot};
+ } else if(vsplinetype.length != 3) abort("vsplinetype must have length 3");
+
+ surface sx=surface(fx,ipt,jpt,usplinetype[0],vsplinetype[0],active);
+ surface sy=surface(fy,ipt,jpt,usplinetype[1],vsplinetype[1],active);
+ surface sz=surface(fz,ipt,jpt,usplinetype[2],vsplinetype[2],active);
+
+ surface s=surface(sx.s.length);
+ for(int k=0; k < sx.s.length; ++k) {
+ triple[][] Q=new triple[4][4];
+ for(int i=0; i < 4 ; ++i) {
+ for(int j=0; j < 4 ; ++j) {
+ Q[i][j]=(sx.s[k].P[i][j].z,sy.s[k].P[i][j].z,sz.s[k].P[i][j].z);
+ }
+ s.s[k]=patch(Q);
+ }
+ }
+ return s;
+}
+
// return the surface described by a real function f over box(a,b),
// interpolated linearly.
surface surface(real f(pair z), pair a, pair b, int nx=nmesh, int ny=nx,
@@ -1741,9 +1817,10 @@ surface surface(real f(pair z), pair a, pair b, int nx=nmesh, int ny=nx,
}
// return the surface described by a real function f over box(a,b),
-// interpolated with splinetype.
+// interpolated with xsplinetype and ysplinetype.
surface surface(real f(pair z), pair a, pair b, int nx=nmesh, int ny=nx,
- splinetype splinetype, bool cond(pair z)=null)
+ splinetype xsplinetype, splinetype ysplinetype=xsplinetype,
+ bool cond(pair z)=null)
{
bool[][] active;
bool all=cond == null;
@@ -1767,7 +1844,7 @@ surface surface(real f(pair z), pair a, pair b, int nx=nmesh, int ny=nx,
if(!all) activei[j]=cond(z);
}
}
- return surface(F,x,y,splinetype,active);
+ return surface(F,x,y,xsplinetype,ysplinetype,active);
}
guide3[][] lift(real f(real x, real y), guide[][] g,
diff --git a/Build/source/utils/asymptote/base/lmfit.asy b/Build/source/utils/asymptote/base/lmfit.asy
new file mode 100644
index 00000000000..d0367547b74
--- /dev/null
+++ b/Build/source/utils/asymptote/base/lmfit.asy
@@ -0,0 +1,838 @@
+/*
+ Copyright (c) 2009 Philipp Stephani
+
+ Permission is hereby granted, free of charge, to any person
+ obtaining a copy of this software and associated documentation files
+ (the "Software"), to deal in the Software without restriction,
+ including without limitation the rights to use, copy, modify, merge,
+ publish, distribute, sublicense, and/or sell copies of the Software,
+ and to permit persons to whom the Software is furnished to do so,
+ subject to the following conditions:
+
+ The above copyright notice and this permission notice shall be
+ included in all copies or substantial portions of the Software.
+
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+ BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+ ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+ CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ SOFTWARE.
+*/
+
+/*
+ This module provides an implementation of the Levenberg--Marquardt
+ (LM) algorithm, converted from the C lmfit routine by Joachim Wuttke
+ (see http://www.messen-und-deuten.de/lmfit/).
+
+ Implementation strategy: Fortunately, Asymptote's syntax is very
+ similar to C, and the original code cleanly separates the
+ customizable parts (user-provided data, output routines, etc.) from
+ the dirty number crunching. Thus, mst of the code was just copied
+ and slightly modified from the original source files. I have
+ amended the lm_data_type structure and the callback routines with a
+ weight array that can be used to provide experimental errors. I
+ have also created two simple wrapper functions.
+*/
+
+
+// copied from the C code
+private real LM_MACHEP = realEpsilon;
+private real LM_DWARF = realMin;
+private real LM_SQRT_DWARF = sqrt(realMin);
+private real LM_SQRT_GIANT = sqrt(realMax);
+private real LM_USERTOL = 30 * LM_MACHEP;
+
+restricted string lm_infmsg[] = {
+ "improper input parameters",
+ "the relative error in the sum of squares is at most tol",
+ "the relative error between x and the solution is at most tol",
+ "both errors are at most tol",
+ "fvec is orthogonal to the columns of the jacobian to machine precision",
+ "number of calls to fcn has reached or exceeded maxcall*(n+1)",
+ "ftol is too small: no further reduction in the sum of squares is possible",
+ "xtol too small: no further improvement in approximate solution x possible",
+ "gtol too small: no further improvement in approximate solution x possible",
+ "not enough memory",
+ "break requested within function evaluation"
+};
+
+restricted string lm_shortmsg[] = {
+ "invalid input",
+ "success (f)",
+ "success (p)",
+ "success (f,p)",
+ "degenerate",
+ "call limit",
+ "failed (f)",
+ "failed (p)",
+ "failed (o)",
+ "no memory",
+ "user break"
+};
+
+
+// copied from the C code and amended with the weight (user_w) array
+struct lm_data_type {
+ real[] user_t;
+ real[] user_y;
+ real[] user_w;
+ real user_func(real user_t_point, real[] par);
+};
+
+
+// Asymptote has no pointer support, so we need reference wrappers for
+// the int and real types
+struct lm_int_type {
+ int val;
+
+ void operator init(int val) {
+ this.val = val;
+ }
+};
+
+
+struct lm_real_type {
+ real val;
+
+ void operator init(real val) {
+ this.val = val;
+ }
+};
+
+
+// copied from the C code; the lm_initialize_control function turned
+// into a constructor
+struct lm_control_type {
+ real ftol;
+ real xtol;
+ real gtol;
+ real epsilon;
+ real stepbound;
+ real fnorm;
+ int maxcall;
+ lm_int_type nfev;
+ lm_int_type info;
+
+ void operator init() {
+ maxcall = 100;
+ epsilon = LM_USERTOL;
+ stepbound = 100;
+ ftol = LM_USERTOL;
+ xtol = LM_USERTOL;
+ gtol = LM_USERTOL;
+ }
+};
+
+
+// copied from the C code
+typedef void lm_evaluate_ftype(real[] par, int m_dat, real[] fvec, lm_data_type data, lm_int_type info);
+typedef void lm_print_ftype(int n_par, real[] par, int m_dat, real[] fvec, lm_data_type data, int iflag, int iter, int nfev);
+
+
+// copied from the C code
+private real SQR(real x) {
+ return x * x;
+}
+
+
+// Asymptote doesn't support pointers to arbitrary array elements, so
+// we provide an offset parameter.
+private real lm_enorm(int n, real[] x, int offset=0) {
+ real s1 = 0;
+ real s2 = 0;
+ real s3 = 0;
+ real x1max = 0;
+ real x3max = 0;
+ real agiant = LM_SQRT_GIANT / n;
+ real xabs, temp;
+
+ for (int i = 0; i < n; ++i) {
+ xabs = fabs(x[offset + i]);
+ if (xabs > LM_SQRT_DWARF && xabs < agiant) {
+ s2 += SQR(xabs);
+ continue;
+ }
+
+ if (xabs > LM_SQRT_DWARF) {
+ if (xabs > x1max) {
+ temp = x1max / xabs;
+ s1 = 1 + s1 * SQR(temp);
+ x1max = xabs;
+ } else {
+ temp = xabs / x1max;
+ s1 += SQR(temp);
+ }
+ continue;
+ }
+ if (xabs > x3max) {
+ temp = x3max / xabs;
+ s3 = 1 + s3 * SQR(temp);
+ x3max = xabs;
+ } else {
+ if (xabs != 0.0) {
+ temp = xabs / x3max;
+ s3 += SQR(temp);
+ }
+ }
+ }
+
+ if (s1 != 0)
+ return x1max * sqrt(s1 + (s2 / x1max) / x1max);
+ if (s2 != 0) {
+ if (s2 >= x3max)
+ return sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));
+ else
+ return sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
+ }
+
+ return x3max * sqrt(s3);
+}
+
+
+// This function calculated the vector whose square sum is to be
+// minimized. We use a slight modification of the original code that
+// includes the weight factor. The user may provide different
+// customizations.
+void lm_evaluate_default(real[] par, int m_dat, real[] fvec, lm_data_type data, lm_int_type info) {
+ for (int i = 0; i < m_dat; ++i) {
+ fvec[i] = data.user_w[i] * (data.user_y[i] - data.user_func(data.user_t[i], par));
+ }
+}
+
+
+// Helper functions to print padded strings and numbers (until
+// Asymptote provides a real printf function)
+private string pad(string str, int count, string pad=" ") {
+ string res = str;
+ while (length(res) < count)
+ res = pad + res;
+ return res;
+}
+
+
+private string pad(int num, int digits, string pad=" ") {
+ return pad(string(num), digits, pad);
+}
+
+
+private string pad(real num, int digits, string pad=" ") {
+ return pad(string(num), digits, pad);
+}
+
+
+// Similar to the C code, also prints weights
+void lm_print_default(int n_par, real[] par, int m_dat, real[] fvec, lm_data_type data, int iflag, int iter, int nfev) {
+ real f, y, t, w;
+ int i;
+
+ if (iflag == 2) {
+ write("trying step in gradient direction");
+ } else if (iflag == 1) {
+ write(format("determining gradient (iteration %d)", iter));
+ } else if (iflag == 0) {
+ write("starting minimization");
+ } else if (iflag == -1) {
+ write(format("terminated after %d evaluations", nfev));
+ }
+
+ write(" par: ", none);
+ for (i = 0; i < n_par; ++i) {
+ write(" " + pad(par[i], 12), none);
+ }
+ write(" => norm: " + pad(lm_enorm(m_dat, fvec), 12));
+
+ if (iflag == -1) {
+ write(" fitting data as follows:");
+ for (i = 0; i < m_dat; ++i) {
+ t = data.user_t[i];
+ y = data.user_y[i];
+ w = data.user_w[i];
+ f = data.user_func(t, par);
+ write(format(" t[%2d]=", i) + pad(t, 12) + " y=" + pad(y, 12) + " w=" + pad(w, 12) + " fit=" + pad(f, 12) + " residue=" + pad(y - f, 12));
+ }
+ }
+}
+
+
+// Prints nothing
+void lm_print_quiet(int n_par, real[] par, int m_dat, real[] fvec, lm_data_type data, int iflag, int iter, int nfev) {
+}
+
+
+// copied from the C code
+private void lm_qrfac(int m, int n, real[] a, bool pivot, int[] ipvt, real[] rdiag, real[] acnorm, real[] wa) {
+ int i, j, k, kmax, minmn;
+ real ajnorm, sum, temp;
+ static real p05 = 0.05;
+
+ for (j = 0; j < n; ++j) {
+ acnorm[j] = lm_enorm(m, a, j * m);
+ rdiag[j] = acnorm[j];
+ wa[j] = rdiag[j];
+ if (pivot)
+ ipvt[j] = j;
+ }
+
+ minmn = min(m, n);
+ for (j = 0; j < minmn; ++j) {
+ while (pivot) {
+ kmax = j;
+ for (k = j + 1; k < n; ++k)
+ if (rdiag[k] > rdiag[kmax])
+ kmax = k;
+ if (kmax == j)
+ break;
+
+ for (i = 0; i < m; ++i) {
+ temp = a[j * m + i];
+ a[j * m + i] = a[kmax * m + i];
+ a[kmax * m + i] = temp;
+ }
+ rdiag[kmax] = rdiag[j];
+ wa[kmax] = wa[j];
+ k = ipvt[j];
+ ipvt[j] = ipvt[kmax];
+ ipvt[kmax] = k;
+
+ break;
+ }
+
+ ajnorm = lm_enorm(m - j, a, j * m + j);
+ if (ajnorm == 0.0) {
+ rdiag[j] = 0;
+ continue;
+ }
+
+ if (a[j * m + j] < 0.0)
+ ajnorm = -ajnorm;
+ for (i = j; i < m; ++i)
+ a[j * m + i] /= ajnorm;
+ a[j * m + j] += 1;
+
+ for (k = j + 1; k < n; ++k) {
+ sum = 0;
+
+ for (i = j; i < m; ++i)
+ sum += a[j * m + i] * a[k * m + i];
+
+ temp = sum / a[j + m * j];
+
+ for (i = j; i < m; ++i)
+ a[k * m + i] -= temp * a[j * m + i];
+
+ if (pivot && rdiag[k] != 0.0) {
+ temp = a[m * k + j] / rdiag[k];
+ temp = max(0.0, 1 - SQR(temp));
+ rdiag[k] *= sqrt(temp);
+ temp = rdiag[k] / wa[k];
+ if (p05 * SQR(temp) <= LM_MACHEP) {
+ rdiag[k] = lm_enorm(m - j - 1, a, m * k + j + 1);
+ wa[k] = rdiag[k];
+ }
+ }
+ }
+
+ rdiag[j] = -ajnorm;
+ }
+}
+
+
+// copied from the C code
+private void lm_qrsolv(int n, real[] r, int ldr, int[] ipvt, real[] diag, real[] qtb, real[] x, real[] sdiag, real[] wa) {
+ static real p25 = 0.25;
+ static real p5 = 0.5;
+
+ int i, kk, j, k, nsing;
+ real qtbpj, sum, temp;
+ real _sin, _cos, _tan, _cot;
+
+ for (j = 0; j < n; ++j) {
+ for (i = j; i < n; ++i)
+ r[j * ldr + i] = r[i * ldr + j];
+ x[j] = r[j * ldr + j];
+ wa[j] = qtb[j];
+ }
+
+ for (j = 0; j < n; ++j) {
+ while (diag[ipvt[j]] != 0.0) {
+ for (k = j; k < n; ++k)
+ sdiag[k] = 0.0;
+ sdiag[j] = diag[ipvt[j]];
+
+ qtbpj = 0.;
+ for (k = j; k < n; ++k) {
+ if (sdiag[k] == 0.)
+ continue;
+ kk = k + ldr * k;
+ if (fabs(r[kk]) < fabs(sdiag[k])) {
+ _cot = r[kk] / sdiag[k];
+ _sin = p5 / sqrt(p25 + p25 * _cot * _cot);
+ _cos = _sin * _cot;
+ } else {
+ _tan = sdiag[k] / r[kk];
+ _cos = p5 / sqrt(p25 + p25 * _tan * _tan);
+ _sin = _cos * _tan;
+ }
+
+ r[kk] = _cos * r[kk] + _sin * sdiag[k];
+ temp = _cos * wa[k] + _sin * qtbpj;
+ qtbpj = -_sin * wa[k] + _cos * qtbpj;
+ wa[k] = temp;
+
+ for (i = k + 1; i < n; ++i) {
+ temp = _cos * r[k * ldr + i] + _sin * sdiag[i];
+ sdiag[i] = -_sin * r[k * ldr + i] + _cos * sdiag[i];
+ r[k * ldr + i] = temp;
+ }
+ }
+ break;
+ }
+
+ sdiag[j] = r[j * ldr + j];
+ r[j * ldr + j] = x[j];
+ }
+
+ nsing = n;
+ for (j = 0; j < n; ++j) {
+ if (sdiag[j] == 0.0 && nsing == n)
+ nsing = j;
+ if (nsing < n)
+ wa[j] = 0;
+ }
+
+ for (j = nsing - 1; j >= 0; --j) {
+ sum = 0;
+ for (i = j + 1; i < nsing; ++i)
+ sum += r[j * ldr + i] * wa[i];
+ wa[j] = (wa[j] - sum) / sdiag[j];
+ }
+
+ for (j = 0; j < n; ++j)
+ x[ipvt[j]] = wa[j];
+}
+
+
+// copied from the C code
+private void lm_lmpar(int n, real[] r, int ldr, int[] ipvt, real[] diag, real[] qtb, real delta, lm_real_type par, real[] x, real[] sdiag, real[] wa1, real[] wa2) {
+ static real p1 = 0.1;
+ static real p001 = 0.001;
+
+ int nsing = n;
+ real parl = 0.0;
+
+ int i, iter, j;
+ real dxnorm, fp, fp_old, gnorm, parc, paru;
+ real sum, temp;
+
+ for (j = 0; j < n; ++j) {
+ wa1[j] = qtb[j];
+ if (r[j * ldr + j] == 0 && nsing == n)
+ nsing = j;
+ if (nsing < n)
+ wa1[j] = 0;
+ }
+ for (j = nsing - 1; j >= 0; --j) {
+ wa1[j] = wa1[j] / r[j + ldr * j];
+ temp = wa1[j];
+ for (i = 0; i < j; ++i)
+ wa1[i] -= r[j * ldr + i] * temp;
+ }
+
+ for (j = 0; j < n; ++j)
+ x[ipvt[j]] = wa1[j];
+
+ iter = 0;
+ for (j = 0; j < n; ++j)
+ wa2[j] = diag[j] * x[j];
+ dxnorm = lm_enorm(n, wa2);
+ fp = dxnorm - delta;
+ if (fp <= p1 * delta) {
+ par.val = 0;
+ return;
+ }
+
+ if (nsing >= n) {
+ for (j = 0; j < n; ++j)
+ wa1[j] = diag[ipvt[j]] * wa2[ipvt[j]] / dxnorm;
+
+ for (j = 0; j < n; ++j) {
+ sum = 0.0;
+ for (i = 0; i < j; ++i)
+ sum += r[j * ldr + i] * wa1[i];
+ wa1[j] = (wa1[j] - sum) / r[j + ldr * j];
+ }
+ temp = lm_enorm(n, wa1);
+ parl = fp / delta / temp / temp;
+ }
+
+ for (j = 0; j < n; ++j) {
+ sum = 0;
+ for (i = 0; i <= j; ++i)
+ sum += r[j * ldr + i] * qtb[i];
+ wa1[j] = sum / diag[ipvt[j]];
+ }
+ gnorm = lm_enorm(n, wa1);
+ paru = gnorm / delta;
+ if (paru == 0.0)
+ paru = LM_DWARF / min(delta, p1);
+
+ par.val = max(par.val, parl);
+ par.val = min(par.val, paru);
+ if (par.val == 0.0)
+ par.val = gnorm / dxnorm;
+
+ for (;; ++iter) {
+ if (par.val == 0.0)
+ par.val = max(LM_DWARF, p001 * paru);
+ temp = sqrt(par.val);
+ for (j = 0; j < n; ++j)
+ wa1[j] = temp * diag[j];
+ lm_qrsolv(n, r, ldr, ipvt, wa1, qtb, x, sdiag, wa2);
+ for (j = 0; j < n; ++j)
+ wa2[j] = diag[j] * x[j];
+ dxnorm = lm_enorm(n, wa2);
+ fp_old = fp;
+ fp = dxnorm - delta;
+
+ if (fabs(fp) <= p1 * delta || (parl == 0.0 && fp <= fp_old && fp_old < 0.0) || iter == 10)
+ break;
+
+ for (j = 0; j < n; ++j)
+ wa1[j] = diag[ipvt[j]] * wa2[ipvt[j]] / dxnorm;
+
+ for (j = 0; j < n; ++j) {
+ wa1[j] = wa1[j] / sdiag[j];
+ for (i = j + 1; i < n; ++i)
+ wa1[i] -= r[j * ldr + i] * wa1[j];
+ }
+ temp = lm_enorm(n, wa1);
+ parc = fp / delta / temp / temp;
+
+ if (fp > 0)
+ parl = max(parl, par.val);
+ else if (fp < 0)
+ paru = min(paru, par.val);
+
+ par.val = max(parl, par.val + parc);
+ }
+}
+
+
+// copied from the C code; the main function
+void lm_lmdif(int m, int n, real[] x, real[] fvec, real ftol, real xtol, real gtol, int maxfev, real epsfcn, real[] diag, int mode, real factor, lm_int_type info, lm_int_type nfev, real[] fjac, int[] ipvt, real[] qtf, real[] wa1, real[] wa2, real[] wa3, real[] wa4, lm_evaluate_ftype evaluate, lm_print_ftype printout, lm_data_type data) {
+ static real p1 = 0.1;
+ static real p5 = 0.5;
+ static real p25 = 0.25;
+ static real p75 = 0.75;
+ static real p0001 = 1.0e-4;
+
+ nfev.val = 0;
+ int iter = 1;
+ lm_real_type par = lm_real_type(0);
+ real delta = 0;
+ real xnorm = 0;
+ real temp = max(epsfcn, LM_MACHEP);
+ real eps = sqrt(temp);
+ int i, j;
+ real actred, dirder, fnorm, fnorm1, gnorm, pnorm, prered, ratio, step, sum, temp1, temp2, temp3;
+
+ if ((n <= 0) || (m < n) || (ftol < 0.0) || (xtol < 0.0) || (gtol < 0.0) || (maxfev <= 0) || (factor <= 0)) {
+ info.val = 0;
+ return;
+ }
+ if (mode == 2) {
+ for (j = 0; j < n; ++j) {
+ if (diag[j] <= 0.0) {
+ info.val = 0;
+ return;
+ }
+ }
+ }
+
+ info.val = 0;
+ evaluate(x, m, fvec, data, info);
+ printout(n, x, m, fvec, data, 0, 0, ++nfev.val);
+ if (info.val < 0)
+ return;
+ fnorm = lm_enorm(m, fvec);
+
+ do {
+ for (j = 0; j < n; ++j) {
+ temp = x[j];
+ step = eps * fabs(temp);
+ if (step == 0.0)
+ step = eps;
+ x[j] = temp + step;
+ info.val = 0;
+ evaluate(x, m, wa4, data, info);
+ printout(n, x, m, wa4, data, 1, iter, ++nfev.val);
+ if (info.val < 0)
+ return;
+ for (i = 0; i < m; ++i)
+ fjac[j * m + i] = (wa4[i] - fvec[i]) / (x[j] - temp);
+ x[j] = temp;
+ }
+
+ lm_qrfac(m, n, fjac, true, ipvt, wa1, wa2, wa3);
+
+ if (iter == 1) {
+ if (mode != 2) {
+ for (j = 0; j < n; ++j) {
+ diag[j] = wa2[j];
+ if (wa2[j] == 0.0)
+ diag[j] = 1.0;
+ }
+ }
+ for (j = 0; j < n; ++j)
+ wa3[j] = diag[j] * x[j];
+ xnorm = lm_enorm(n, wa3);
+ delta = factor * xnorm;
+ if (delta == 0.0)
+ delta = factor;
+ }
+
+ for (i = 0; i < m; ++i)
+ wa4[i] = fvec[i];
+
+ for (j = 0; j < n; ++j) {
+ temp3 = fjac[j * m + j];
+ if (temp3 != 0.0) {
+ sum = 0;
+ for (i = j; i < m; ++i)
+ sum += fjac[j * m + i] * wa4[i];
+ temp = -sum / temp3;
+ for (i = j; i < m; ++i)
+ wa4[i] += fjac[j * m + i] * temp;
+ }
+ fjac[j * m + j] = wa1[j];
+ qtf[j] = wa4[j];
+ }
+
+ gnorm = 0;
+ if (fnorm != 0) {
+ for (j = 0; j < n; ++j) {
+ if (wa2[ipvt[j]] == 0) continue;
+ sum = 0.0;
+ for (i = 0; i <= j; ++i)
+ sum += fjac[j * m + i] * qtf[i] / fnorm;
+ gnorm = max(gnorm, fabs(sum / wa2[ipvt[j]]));
+ }
+ }
+
+ if (gnorm <= gtol) {
+ info.val = 4;
+ return;
+ }
+
+ if (mode != 2) {
+ for (j = 0; j < n; ++j)
+ diag[j] = max(diag[j], wa2[j]);
+ }
+
+ do {
+ lm_lmpar(n, fjac, m, ipvt, diag, qtf, delta, par, wa1, wa2, wa3, wa4);
+
+ for (j = 0; j < n; ++j) {
+ wa1[j] = -wa1[j];
+ wa2[j] = x[j] + wa1[j];
+ wa3[j] = diag[j] * wa1[j];
+ }
+ pnorm = lm_enorm(n, wa3);
+
+ if (nfev.val <= 1 + n)
+ delta = min(delta, pnorm);
+
+ info.val = 0;
+ evaluate(wa2, m, wa4, data, info);
+ printout(n, x, m, wa4, data, 2, iter, ++nfev.val);
+ if (info.val < 0)
+ return;
+
+ fnorm1 = lm_enorm(m, wa4);
+
+ if (p1 * fnorm1 < fnorm)
+ actred = 1 - SQR(fnorm1 / fnorm);
+ else
+ actred = -1;
+
+ for (j = 0; j < n; ++j) {
+ wa3[j] = 0;
+ for (i = 0; i <= j; ++i)
+ wa3[i] += fjac[j * m + i] * wa1[ipvt[j]];
+ }
+ temp1 = lm_enorm(n, wa3) / fnorm;
+ temp2 = sqrt(par.val) * pnorm / fnorm;
+ prered = SQR(temp1) + 2 * SQR(temp2);
+ dirder = -(SQR(temp1) + SQR(temp2));
+
+ ratio = prered != 0 ? actred / prered : 0;
+
+ if (ratio <= p25) {
+ if (actred >= 0.0)
+ temp = p5;
+ else
+ temp = p5 * dirder / (dirder + p5 * actred);
+ if (p1 * fnorm1 >= fnorm || temp < p1)
+ temp = p1;
+ delta = temp * min(delta, pnorm / p1);
+ par.val /= temp;
+ } else if (par.val == 0.0 || ratio >= p75) {
+ delta = pnorm / p5;
+ par.val *= p5;
+ }
+
+ if (ratio >= p0001) {
+ for (j = 0; j < n; ++j) {
+ x[j] = wa2[j];
+ wa2[j] = diag[j] * x[j];
+ }
+ for (i = 0; i < m; ++i)
+ fvec[i] = wa4[i];
+ xnorm = lm_enorm(n, wa2);
+ fnorm = fnorm1;
+ ++iter;
+ }
+
+ info.val = 0;
+ if (fabs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1)
+ info.val = 1;
+ if (delta <= xtol * xnorm)
+ info.val += 2;
+ if (info.val != 0)
+ return;
+
+ if (nfev.val >= maxfev)
+ info.val = 5;
+ if (fabs(actred) <= LM_MACHEP && prered <= LM_MACHEP && p5 * ratio <= 1)
+ info.val = 6;
+ if (delta <= LM_MACHEP * xnorm)
+ info.val = 7;
+ if (gnorm <= LM_MACHEP)
+ info.val = 8;
+ if (info.val != 0)
+ return;
+ } while (ratio < p0001);
+ } while (true);
+}
+
+
+// copied from the C code; wrapper of lm_lmdif
+void lm_minimize(int m_dat, int n_par, real[] par, lm_evaluate_ftype evaluate, lm_print_ftype printout, lm_data_type data, lm_control_type control) {
+ int n = n_par;
+ int m = m_dat;
+
+ real[] fvec = new real[m];
+ real[] diag = new real[n];
+ real[] qtf = new real[n];
+ real[] fjac = new real[n * m];
+ real[] wa1 = new real[n];
+ real[] wa2 = new real[n];
+ real[] wa3 = new real[n];
+ real[] wa4 = new real[m];
+ int[] ipvt = new int[n];
+
+ control.info.val = 0;
+ control.nfev.val = 0;
+
+ lm_lmdif(m, n, par, fvec, control.ftol, control.xtol, control.gtol, control.maxcall * (n + 1), control.epsilon, diag, 1, control.stepbound, control.info, control.nfev, fjac, ipvt, qtf, wa1, wa2, wa3, wa4, evaluate, printout, data);
+
+ printout(n, par, m, fvec, data, -1, 0, control.nfev.val);
+ control.fnorm = lm_enorm(m, fvec);
+ if (control.info.val < 0)
+ control.info.val = 10;
+}
+
+
+// convenience functions; wrappers of lm_minimize
+struct FitControl {
+ real squareSumTolerance;
+ real approximationTolerance;
+ real desiredOrthogonality;
+ real epsilon;
+ real stepBound;
+ int maxIterations;
+ bool verbose;
+
+ void operator init(real squareSumTolerance, real approximationTolerance, real desiredOrthogonality, real epsilon, real stepBound, int maxIterations, bool verbose) {
+ this.squareSumTolerance = squareSumTolerance;
+ this.approximationTolerance = approximationTolerance;
+ this.desiredOrthogonality = desiredOrthogonality;
+ this.epsilon = epsilon;
+ this.stepBound = stepBound;
+ this.maxIterations = maxIterations;
+ this.verbose = verbose;
+ }
+
+ FitControl copy() {
+ FitControl result = new FitControl;
+ result.squareSumTolerance = this.squareSumTolerance;
+ result.approximationTolerance = this.approximationTolerance;
+ result.desiredOrthogonality = this.desiredOrthogonality;
+ result.epsilon = this.epsilon;
+ result.stepBound = this.stepBound;
+ result.maxIterations = this.maxIterations;
+ result.verbose = this.verbose;
+ return result;
+ }
+};
+
+FitControl defaultControl;
+defaultControl.squareSumTolerance = LM_USERTOL;
+defaultControl.approximationTolerance = LM_USERTOL;
+defaultControl.desiredOrthogonality = LM_USERTOL;
+defaultControl.epsilon = LM_USERTOL;
+defaultControl.stepBound = 100;
+defaultControl.maxIterations = 100;
+defaultControl.verbose = false;
+
+
+struct FitResult {
+ restricted real norm;
+ restricted int iterations;
+ restricted int status;
+
+ void operator init(real norm, int status, int iterations) {
+ this.norm = norm;
+ this.status = status;
+ this.iterations = iterations;
+ }
+};
+
+
+// Fits data points (xdata, ydata ± errors) to the given function using the given parameters.
+FitResult fit(real[] xdata, real[] ydata, real[] errors, real function(real[], real), real[] parameters, FitControl control=defaultControl) {
+ int m_dat = min(xdata.length, ydata.length);
+ int n_par = parameters.length;
+ lm_evaluate_ftype evaluate = lm_evaluate_default;
+ lm_print_ftype printout = control.verbose ? lm_print_default : lm_print_quiet;
+
+ lm_data_type data;
+ data.user_t = xdata;
+ data.user_y = ydata;
+ data.user_w = 1 / errors;
+ data.user_func = new real(real x, real[] params) {
+ return function(params, x);
+ };
+
+ lm_control_type ctrl;
+ ctrl.ftol = control.squareSumTolerance;
+ ctrl.xtol = control.approximationTolerance;
+ ctrl.gtol = control.desiredOrthogonality;
+ ctrl.epsilon = control.epsilon;
+ ctrl.stepbound = control.stepBound;
+ ctrl.maxcall = control.maxIterations;
+
+ lm_minimize(m_dat, n_par, parameters, evaluate, printout, data, ctrl);
+
+ return FitResult(ctrl.fnorm, ctrl.nfev.val, ctrl.info.val);
+}
+
+
+// Fits data points (xdata, ydata) to the given function using the given parameters.
+FitResult fit(real[] xdata, real[] ydata, real function(real[], real), real[] parameters, FitControl control=defaultControl) {
+ return fit(xdata, ydata, array(min(xdata.length, ydata.length), 1.0), function, parameters, control);
+}
+
diff --git a/Build/source/utils/asymptote/base/math.asy b/Build/source/utils/asymptote/base/math.asy
index f204f746349..85d41d86d07 100644
--- a/Build/source/utils/asymptote/base/math.asy
+++ b/Build/source/utils/asymptote/base/math.asy
@@ -148,7 +148,7 @@ real[][] zero(int n, int m)
return M;
}
-real[][] operator + (real[][] a, real[][] b)
+real[][] operator +(real[][] a, real[][] b)
{
int n=a.length;
real[][] m=new real[n][];
@@ -157,7 +157,7 @@ real[][] operator + (real[][] a, real[][] b)
return m;
}
-real[][] operator - (real[][] a, real[][] b)
+real[][] operator -(real[][] a, real[][] b)
{
int n=a.length;
real[][] m=new real[n][];
@@ -166,26 +166,7 @@ real[][] operator - (real[][] a, real[][] b)
return m;
}
-private string incommensurate=
- "Multiplication of incommensurate matrices is undefined";
-
-real[] operator * (real[] b, real[][] a)
-{
- int nb=b.length;
- if(nb != a.length)
- abort(incommensurate);
- int na0=a[0].length;
- real[] m=new real[na0];
- for(int j=0; j < na0; ++j) {
- real sum;
- for(int k=0; k < nb; ++k)
- sum += b[k]*a[k][j];
- m[j]=sum;
- }
- return m;
-}
-
-real[][] operator * (real[][] a, real b)
+real[][] operator *(real[][] a, real b)
{
int n=a.length;
real[][] m=new real[n][];
@@ -194,12 +175,12 @@ real[][] operator * (real[][] a, real b)
return m;
}
-real[][] operator * (real b, real[][] a)
+real[][] operator *(real b, real[][] a)
{
return a*b;
}
-real[][] operator / (real[][] a, real b)
+real[][] operator /(real[][] a, real b)
{
return a*(1/b);
}
diff --git a/Build/source/utils/asymptote/base/ode.asy b/Build/source/utils/asymptote/base/ode.asy
index 5135dd9443f..643c56feec7 100644
--- a/Build/source/utils/asymptote/base/ode.asy
+++ b/Build/source/utils/asymptote/base/ode.asy
@@ -1,15 +1,260 @@
-real euler(real y, real f(real x, real y), real a, real b=a, int n=0,
- real h=0, bool dynamic=false, real tolmin=0, real tolmax=0)
+real stepfactor=2.0; // Maximum dynamic step size adjustment factor.
+
+struct RKTableau
+{
+ int order;
+ real[] steps;
+ real[][] weights;
+ real[] highOrderWeights;
+ real[] lowOrderWeights;
+ real pgrow;
+ real pshrink;
+
+ void operator init(int order, real[][] weights, real[] highOrderWeights,
+ real[] lowOrderWeights=new real[],
+ real[] steps=sequence(new real(int i) {
+ return sum(weights[i]);},weights.length)) {
+ this.order=order;
+ this.steps=steps;
+ this.weights=weights;
+ this.highOrderWeights=highOrderWeights;
+ this.lowOrderWeights=lowOrderWeights;
+ pgrow=(order > 0) ? 1/order : 0;
+ pshrink=(order > 1) ? 1/(order-1) : pgrow;
+ }
+}
+
+// First-Order Euler
+RKTableau Euler=RKTableau(1,new real[][],
+ new real[] {1});
+
+// Second-Order Runge-Kutta
+RKTableau RK2=RKTableau(2,new real[][] {{1/2}},
+ new real[] {0,1});
+
+// Second-Order Predictor-Corrector
+RKTableau PC=RKTableau(2,new real[][] {{1}},
+ new real[] {1/2,1/2});
+
+// Third-Order Classical Runge-Kutta
+RKTableau RK3=RKTableau(3,new real[][] {{1/2},{-1,2}},
+ new real[] {1/6,2/3,1/6});
+
+// Third-Order Bogacki-Shampine Runge-Kutta
+RKTableau RK3BS=RKTableau(3,new real[][] {{1/2},{0,3/4}},
+ new real[] {2/9,1/3,4/9}, // 3rd order
+ new real[] {7/24,1/4,1/3,1/8}); // 2nd order
+
+// Fourth-Order Classical Runge-Kutta
+RKTableau RK4=RKTableau(4,new real[][] {{1/2},{0,1/2},{0,0,1}},
+ new real[] {1/6,1/3,1/3,1/6});
+
+// Fifth-Order Cash-Karp Runge-Kutta
+RKTableau RK5CK=RKTableau(5,new real[][] {{1/5},
+ {3/40,9/40},
+ {3/10,-9/10,6/5},
+ {-11/54,5/2,-70/27,35/27},
+ {1631/55296,175/512,575/13824,
+ 44275/110592,253/4096}},
+ new real[] {37/378,0,250/621,125/594,
+ 0,512/1771}, // 5th order
+ new real[] {2825/27648,0,18575/48384,13525/55296,
+ 277/14336,1/4}); // 4th order
+
+// Fifth-Order Fehlberg Runge-Kutta
+RKTableau RK5F=RKTableau(5,new real[][] {{1/4},
+ {3/32,9/32},
+ {1932/2197,-7200/2197,7296/2197},
+ {439/216,-8,3680/513,-845/4104},
+ {-8/27,2,-3544/2565,1859/4104,
+ -11/40}},
+ new real[] {16/135,0,6656/12825,28561/56430,-9/50,2/55}, // 5th order
+ new real[] {25/216,0,1408/2565,2197/4104,-1/5,0}); // 4th order
+
+// Fifth-Order Dormand-Prince Runge-Kutta
+RKTableau RK5DP=RKTableau(5,new real[][] {{1/5},
+ {3/40,9/40},
+ {44/45,-56/15,32/9},
+ {19372/6561,-25360/2187,64448/6561,
+ -212/729},
+ {9017/3168,-355/33,46732/5247,49/176,
+ -5103/18656}},
+ new real[] {35/384,0,500/1113,125/192,-2187/6784,
+ 11/84}, // 5th order
+ new real[] {5179/57600,0,7571/16695,393/640,
+ -92097/339200,187/2100,1/40}); // 4th order
+
+real error(real error, real initial, real norm, real lowOrder, real diff)
+{
+ if(initial != 0.0 && lowOrder != initial) {
+ static real epsilon=realMin/realEpsilon;
+ real denom=max(abs(norm),abs(initial))+epsilon;
+ return max(error,max(abs(diff)/denom));
+ }
+ return error;
+}
+
+real adjust(real h, real error, real t, real tolmin, real tolmax,
+ real dtmin, real dtmax, RKTableau tableau, bool verbose=true)
+{
+ real dt=h;
+ void report(real t) {
+ if(h != dt)
+ write("Time step changed from "+(string) dt+" to "+(string) h+" at t="+
+ (string) t+".");
+ }
+ if(error > tolmax) {
+ h=max(h*max((tolmin/error)^tableau.pshrink,1/stepfactor),dtmin);
+ if(verbose) report(t);
+ return h;
+ }
+ if(error > 0 && error < tolmin) {
+ h=min(h*min((tolmin/error)^tableau.pgrow,stepfactor),dtmax);
+ if(verbose) report(t+dt);
+ }
+ return h;
+}
+
+// Integrate dy/dt=f(t,y) from a to b using initial conditions y,
+// specifying either the step size h or the number of steps n.
+real integrate(real y, real f(real t, real y), real a, real b=a, real h=0,
+ int n=0, bool dynamic=false, real tolmin=0, real tolmax=0,
+ real dtmin=0, real dtmax=realMax, RKTableau tableau,
+ bool verbose=false)
+{
+ if(h == 0) {
+ if(b == a) return y;
+ if(n == 0) abort("Either n or h must be specified");
+ else h=(b-a)/n;
+ }
+ real t=a;
+ real f0;
+ bool fsal=tableau.lowOrderWeights.length > tableau.highOrderWeights.length;
+ if(fsal) f0=f(t,y);
+ if(tableau.lowOrderWeights.length == 0) dynamic=false;
+
+ while(t < b) {
+ real[] predictions={fsal ? f0 : f(t,y)};
+ for(int i=0; i < tableau.steps.length; ++i)
+ predictions.push(f(t+h*tableau.steps[i],
+ y+h*dot(tableau.weights[i],predictions)));
+
+ real highOrder=h*dot(tableau.highOrderWeights,predictions);
+ if(dynamic) {
+ real f1;
+ if(fsal) {
+ f1=f(t+h,y+highOrder);
+ predictions.push(f1);
+ }
+ real lowOrder=h*dot(tableau.lowOrderWeights,predictions);
+ real error;
+ error=error(error,y,y+highOrder,y+lowOrder,highOrder-lowOrder);
+ real dt=h;
+ h=adjust(h,error,t,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose);
+ if(h >= dt) {
+ t += dt;
+ y += highOrder;
+ f0=f1;
+ }
+ } else {
+ t += h;
+ y += highOrder;
+ }
+ h=min(h,b-t);
+ if(t >= b || t+h == t) break;
+ }
+ return y;
+}
+
+// Integrate a set of equations, dy/dt=f(t,y), from a to b using initial
+// conditions y, specifying either the step size h or the number of steps n.
+real[] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a,
+ real h=0, int n=0, bool dynamic=false,
+ real tolmin=0, real tolmax=0, real dtmin=0, real dtmax=realMax,
+ RKTableau tableau, bool verbose=false)
{
if(h == 0) {
if(b == a) return y;
if(n == 0) abort("Either n or h must be specified");
else h=(b-a)/n;
}
- real x=a;
- for(int i=0; i < n; ++i) {
- y += h*f(x,y);
- x += h;
+ real[] y=copy(y);
+ real t=a;
+ real[] f0;
+ bool fsal=tableau.lowOrderWeights.length > tableau.highOrderWeights.length;
+ if(fsal) f0=f(t,y);
+ if(tableau.lowOrderWeights.length == 0) dynamic=false;
+
+ while(t < b) {
+ real[][] predictions={fsal ? f0 : f(t,y)};
+ for(int i=0; i < tableau.steps.length; ++i)
+ predictions.push(f(t+h*tableau.steps[i],
+ y+h*tableau.weights[i]*predictions));
+
+ real[] highOrder=h*tableau.highOrderWeights*predictions;
+ if(dynamic) {
+ real[] f1;
+ if(fsal) {
+ f1=f(t+h,y+highOrder);
+ predictions.push(f1);
+ }
+ real[] lowOrder=h*tableau.lowOrderWeights*predictions;
+ real error;
+ for(int i=0; i < y.length; ++i)
+ error=error(error,y[i],y[i]+highOrder[i],y[i]+lowOrder[i],
+ highOrder[i]-lowOrder[i]);
+ real dt=h;
+ h=adjust(h,error,t,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose);
+ if(h >= dt) {
+ t += dt;
+ y += highOrder;
+ f0=f1;
+ }
+ } else {
+ t += h;
+ y += highOrder;
+ }
+ h=min(h,b-t);
+ if(t >= b || t+h == t) break;
}
return y;
}
+
+real[][] finiteDifferenceJacobian(real[] f(real[]), real[] t,
+ real[] h=sqrtEpsilon*abs(t))
+{
+ real[] ft=f(t);
+ real[][] J=new real[t.length][ft.length];
+ real[] ti=copy(t);
+ real tlast=ti[0];
+ ti[0] += h[0];
+ J[0]=(f(ti)-ft)/h[0];
+ for(int i=1; i < t.length; ++i) {
+ ti[i-1]=tlast;
+ tlast=ti[i];
+ ti[i] += h[i];
+ J[i]=(f(ti)-ft)/h[i];
+ }
+ return transpose(J);
+}
+
+// Solve simultaneous nonlinear system by Newton's method.
+real[] newton(int iterations=100, real[] f(real[]), real[][] jacobian(real[]),
+ real[] t)
+{
+ real[] t=copy(t);
+ for(int i=0; i < iterations; ++i)
+ t += solve(jacobian(t),-f(t));
+ return t;
+}
+
+real[] solveBVP(real[] f(real, real[]), real a, real b=a, real h=0, int n=0,
+ real[] initial(real[]), real[] discrepancy(real[]),
+ real[] guess, RKTableau tableau, int iterations=100)
+{
+ real[] g(real[] t) {
+ return discrepancy(integrate(initial(t),f,a,b,h,n,tableau));
+ }
+ real[][] jacobian(real[] t) {return finiteDifferenceJacobian(g,t);}
+ return initial(newton(iterations,g,jacobian,guess));
+}
diff --git a/Build/source/utils/asymptote/base/plain.asy b/Build/source/utils/asymptote/base/plain.asy
index a918a6bdc09..8424c767ffa 100644
--- a/Build/source/utils/asymptote/base/plain.asy
+++ b/Build/source/utils/asymptote/base/plain.asy
@@ -38,6 +38,8 @@ include plain_markers;
include plain_arrows;
include plain_debugger;
+pair viewportmargin=(0,0); // Horizontal and vertical 3D viewport margins.
+
typedef void exitfcn();
bool needshipout() {
diff --git a/Build/source/utils/asymptote/base/plain_strings.asy b/Build/source/utils/asymptote/base/plain_strings.asy
index e3d0a797d86..0fe59ba5a63 100644
--- a/Build/source/utils/asymptote/base/plain_strings.asy
+++ b/Build/source/utils/asymptote/base/plain_strings.asy
@@ -181,9 +181,9 @@ string math(real x)
return math((string) x);
}
-string format(real x)
+string format(real x, string locale="")
{
- return format(defaultformat,x);
+ return format(defaultformat,x,locale);
}
string phantom(string s)
diff --git a/Build/source/utils/asymptote/base/slide.asy b/Build/source/utils/asymptote/base/slide.asy
index 5b3596d02e6..b154c1d4746 100644
--- a/Build/source/utils/asymptote/base/slide.asy
+++ b/Build/source/utils/asymptote/base/slide.asy
@@ -519,7 +519,6 @@ void bibliographystyle(string name)
{
settings.twice=true;
settings.keepaux=true;
- delete(outprefix()+"_.aux");
texpreamble("\bibliographystyle{"+name+"}");
}
diff --git a/Build/source/utils/asymptote/base/splinetype.asy b/Build/source/utils/asymptote/base/splinetype.asy
index 366f91ac370..08263ddf4ca 100644
--- a/Build/source/utils/asymptote/base/splinetype.asy
+++ b/Build/source/utils/asymptote/base/splinetype.asy
@@ -2,6 +2,7 @@ typedef real[] splinetype(real[], real[]);
restricted real[] defaultspline(real[] x, real[] y);
restricted real[] Spline(real[] x, real[] y);
+restricted splinetype[] Spline;
string morepoints="interpolation requires at least 2 points";
string differentlengths="arrays have different lengths";
@@ -63,7 +64,8 @@ real[] periodic(real[] x, real[] y)
{
int n=x.length;
checklengths(n,y.length);
- if(y[n-1] != y[0]) abort("function values are not periodic");
+ if(abs(y[n-1]-y[0]) > sqrtEpsilon*max(abs(y)))
+ abort("function values are not periodic");
real[] d;
if(n > 2) {
real[] a=new real[n-1];
diff --git a/Build/source/utils/asymptote/base/three.asy b/Build/source/utils/asymptote/base/three.asy
index 0af0d2c3bdf..d454d032b23 100644
--- a/Build/source/utils/asymptote/base/three.asy
+++ b/Build/source/utils/asymptote/base/three.asy
@@ -16,11 +16,10 @@ real defaultgranularity=0;
real linegranularity=0.01;
real tubegranularity=0.003;
real dotgranularity=0.0001;
-pair viewportmargin=0; // Horizontal and vertical viewport margins.
real viewportfactor=1.002; // Factor used to expand orthographic viewport.
-real viewportpadding=1.2; // Offset used to expand PRC viewport.
+real viewportpadding=1; // Offset used to expand PRC viewport.
real angleprecision=1e-3; // Precision for centering perspective projections.
-real anglefactor=max(1.005,1+angleprecision);
+real anglefactor=max(1.01,1+angleprecision);
// Factor used to expand perspective viewport.
string defaultembed3Doptions;
@@ -2178,7 +2177,7 @@ projection perspective(string s)
return P;
}
-private string format(real x)
+private string Format(real x)
{
// Work around movie15.sty division by zero bug;
// e.g. u=unit((1e-10,1e-10,0.9));
@@ -2187,14 +2186,14 @@ private string format(real x)
return format("%.18f",x,"C");
}
-private string format(triple v, string sep=" ")
+private string Format(triple v, string sep=" ")
{
- return format(v.x)+sep+format(v.y)+sep+format(v.z);
+ return Format(v.x)+sep+Format(v.y)+sep+Format(v.z);
}
-private string format(real[] c)
+private string Format(real[] c)
{
- return format((c[0],c[1],c[2]));
+ return Format((c[0],c[1],c[2]));
}
private string[] file3;
@@ -2227,8 +2226,8 @@ string lightscript(light light) {
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';
+ 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 +"
@@ -2255,14 +2254,9 @@ void writeJavaScript(string name, string preamble, string script)
file3.push(name);
}
-pair viewportmargin(real width, real height)
+pair viewportmargin(pair lambda)
{
- pair viewportmargin=viewportmargin;
- real xmargin=viewportmargin.x;
- real ymargin=viewportmargin.y;
- if(xmargin <= 0) xmargin=max(0.5*(viewportsize.x-width),0);
- if(ymargin <= 0) ymargin=max(0.5*(viewportsize.y-height),0);
- return (xmargin,ymargin);
+ return maxbound(0.5*(viewportsize-lambda),viewportmargin);
}
string embed3D(string label="", string text=label, string prefix,
@@ -2284,7 +2278,7 @@ string embed3D(string label="", string text=label, string prefix,
real viewplanesize;
if(P.infinity) {
triple lambda=max3(f)-min3(f);
- pair margin=viewportpadding*(1,1)+viewportmargin(lambda.x,lambda.y);
+ pair margin=viewportpadding*(1,1)+viewportmargin((lambda.x,lambda.y));
viewplanesize=(max(lambda.x+2*margin.x,lambda.y+2*margin.y))/cm;
} else
if(!P.absolute) angle=2*aTan(Tan(0.5*angle)-viewportpadding/P.target.z);
@@ -2314,12 +2308,12 @@ string embed3D(string label="", string text=label, string prefix,
options3 += ",poster";
options3 += ",text={"+text+"},label="+label+
",toolbar="+(settings.toolbar ? "true" : "false")+
- ",3Daac="+format(P.absolute ? P.angle : angle)+
- ",3Dc2c="+format(u)+
- ",3Dcoo="+format(P.target/cm)+
- ",3Droll="+format(roll)+
- ",3Droo="+format(abs(v))+
- ",3Dbg="+format(light.background());
+ ",3Daac="+Format(P.absolute ? P.angle : angle)+
+ ",3Dc2c="+Format(u)+
+ ",3Dcoo="+Format(P.target/cm)+
+ ",3Droll="+Format(roll)+
+ ",3Droo="+Format(abs(v))+
+ ",3Dbg="+Format(light.background());
if(options != "") options3 += ","+options;
if(name != "") options3 += ",3Djscript="+stripdirectory(name);
@@ -2406,7 +2400,7 @@ object embed(string label="", string text=label,
pair m2=pic2.min(s);
pair M2=pic2.max(s);
pair lambda=M2-m2;
- pair viewportmargin=viewportmargin(lambda.x,lambda.y);
+ pair viewportmargin=viewportmargin(lambda);
real width=ceil(lambda.x+2*viewportmargin.x);
real height=ceil(lambda.y+2*viewportmargin.y);
@@ -2440,7 +2434,7 @@ object embed(string label="", string text=label,
triple m=min3(f);
triple M=max3(f);
triple lambda=M-m;
- viewportmargin=viewportmargin(lambda.x,lambda.y);
+ viewportmargin=viewportmargin((lambda.x,lambda.y));
width=lambda.x+2*viewportmargin.x;
height=lambda.y+2*viewportmargin.y;
@@ -2542,7 +2536,6 @@ object embed(string label="", string text=label,
if(P.infinity) {
triple margin=(viewportfactor-1.0)*(abs(M.x-m.x),abs(M.y-m.y),0)
+(viewportmargin.x,viewportmargin.y,0);
-
M += margin;
m -= margin;
} else if(M.z >= 0) abort("camera too close");
diff --git a/Build/source/utils/asymptote/base/three_surface.asy b/Build/source/utils/asymptote/base/three_surface.asy
index 3c49b382133..597104ad202 100644
--- a/Build/source/utils/asymptote/base/three_surface.asy
+++ b/Build/source/utils/asymptote/base/three_surface.asy
@@ -1191,11 +1191,11 @@ void label(picture pic=currentpicture, Label L, triple position,
pic.add(new void(frame f, transform3 t, picture pic, projection P) {
// Handle relative projected 3D alignments.
Label L=L.copy();
+ triple v=t*position;
if(!align.is3D && L.align.relative && L.align.dir3 != O &&
determinant(P.t) != 0)
- L.align(L.align.dir*unit(project(L.align.dir3,P.t)));
+ L.align(L.align.dir*unit(project(v+L.align.dir3,P.t)-project(v,P.t)));
- triple v=t*position;
if(L.defaulttransform3)
L.T3=transform3(P);
if(is3D())
@@ -1234,7 +1234,7 @@ void label(picture pic=currentpicture, Label L, path3 g, align align=NoAlign,
if(L.align.default) {
align a;
a.init(-I*(position <= sqrtEpsilon ? S :
- position >= length(g)-sqrtEpsilon ? N : E),relative=true);
+ position >= length(g)-sqrtEpsilon ? N : E),relative=true);
a.dir3=dir(g,position); // Pass 3D direction via unused field.
L.align(a);
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-18 01:24:18 +0000