diff options
author | Karl Berry <karl@freefriends.org> | 2013-04-07 18:19:31 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2013-04-07 18:19:31 +0000 |
commit | 752012c605d34cd943795527a9738475a6958fcc (patch) | |
tree | 4ee06acdd8333a662c2d6f6ef716235053468f55 /Master/texmf/asymptote/lmfit.asy | |
parent | 9789d09132f18a838e84f041b4b3aff28d3426ec (diff) |
texmf -> texmf-dist: start with unique dirs from texmf
git-svn-id: svn://tug.org/texlive/trunk@29712 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf/asymptote/lmfit.asy')
-rw-r--r-- | Master/texmf/asymptote/lmfit.asy | 881 |
1 files changed, 0 insertions, 881 deletions
diff --git a/Master/texmf/asymptote/lmfit.asy b/Master/texmf/asymptote/lmfit.asy deleted file mode 100644 index e70421074e5..00000000000 --- a/Master/texmf/asymptote/lmfit.asy +++ /dev/null @@ -1,881 +0,0 @@ -/* - 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. -*/ - -/* - Fitting $n$ data points $(x_1, y_1 \pm \Delta y_1), \dots, (x_n, y_n \pm \Delta y_n)$ - to a function $f$ that depends on $m$ parameters $a_1, \dots, a_m$ means minimizing - the least-squares sum - % - \begin{equation*} - \sum_{i = 1}^n \left( \frac{y_i - f(a_1, \dots, a_m; x_i)}{\Delta y_i} \right)^2 - \end{equation*} - % - with respect to the parameters $a_1, \dots, a_m$. -*/ - -/* - 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 - -/* - The structure FitControl specifies various control parameters. -*/ -struct FitControl { - real squareSumTolerance; // relative error desired in the sum of squares - real approximationTolerance; // relative error between last two approximations - real desiredOrthogonality; // orthogonality desired between the residue vector and its derivatives - real epsilon; // step used to calculate the jacobian - real stepBound; // initial bound to steps in the outer loop - int maxIterations; // maximum number of iterations - bool verbose; // whether to print detailed information about every iteration, or nothing - - void operator init(real squareSumTolerance=LM_USERTOL, real approximationTolerance=LM_USERTOL, real desiredOrthogonality=LM_USERTOL, real epsilon=LM_USERTOL, real stepBound=100, int maxIterations=100, bool verbose=false) { - 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 operator init() { - return FitControl(); -} - -FitControl defaultControl; - - -/* - Upon returning, this structure provides information about the fit. -*/ -struct FitResult { - real norm; // norm of the residue vector - int iterations; // actual number of iterations - int status; // status of minimization - - void operator init(real norm, int iterations, int status) { - this.norm = norm; - this.iterations = iterations; - this.status = status; - } -}; - - -/* - Fits data points to a function that depends on some parameters. - - Parameters: - - xdata: Array of x values. - - ydata: Array of y values. - - errors: Array of experimental errors; each element must be strictly positive - - function: Fit function. - - parameters: Parameter array. Before calling fit(), this must contain the initial guesses for the parameters. - Upon return, it will contain the solution parameters. - - control: object of type FitControl that controls various aspects of the fitting procedure. - - Returns: - An object of type FitResult that conveys information about the fitting process. -*/ -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 to a function that depends on some parameters. - - Parameters: - - xdata: Array of x values. - - ydata: Array of y values. - - function: Fit function. - - parameters: Parameter array. Before calling fit(), this must contain the initial guesses for the parameters. - Upon return, it will contain the solution parameters. - - control: object of type FitControl that controls various aspects of the fitting procedure. - - Returns: - An object of type FitResult that conveys information about the fitting process. -*/ -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); -} - |