\ \ curvefit.4th \ \ Nonlinear least-squares curve fitting routine from P.R. Bevington, \ "Data Reduction and Error Analysis for the Physical Sciences", \ translated to kForth by K. Myneni, 9-21-1998 \ \ kForth requires the following files: \ \ ans-words.4th \ fsl-util.4th \ dynmem.4th \ gaussj.4th \ \ Usage: \ \ First include the source file for the array and matrix words \ (fsl-util.4th), then define a fitting function which has the \ following stack diagram: \ \ ( rx a -- ry ) or ( a -- ) ( F: rx -- ry ) \ \ where rx is the x value at which the function is to be \ evaluated and 'a' is the address of the array containing \ the parameter values. The computed value for y is \ returned on the stack (or fpstack). \ \ Input data (x, y) is passed to the fitter in two FSL arrays, \ one containing the x values, the other containing the y values. \ Starting parameters (best guess) for the fit, and reasonable \ increments for the search, are also passed to the fitter. \ A function parameter can be held fixed by setting its increment \ to zero. \ \ Revised: \ \ 2003-12-17 fix problem with using zero parameter increment KM \ 2004-12-02 increased max curvature matrix size to 24x24 KM \ 2006-03-17 modified original version to use FSL-style arrays \ and matrices KM \ 2007-09-21 this file has been renamed from curvefit-fsl.4th \ to curvefit.4th; the original curvefit.4th is now \ obsolete, since we have deprecated the matrix package, \ matrix.4th KM \ 2009-12-15 added error check in CURFIT to ensure that the \ number of terms passed to the routine does not exceed \ MAX_PARAMETERS km \ 2023-11-14 Extensive revisions: set parameter uncertainty \ to zero if parameter increment is zero; moved \ initialization stuff from CURFIT to INIT-CURVEFIT; \ FLAMDA is renamed to LAMBDA, and initialized to 0.001 \ at the start of the fit by INIT-CURVEFIT; add default \ tolerance for fractional change in chi-square for \ CURFIT termination and SET-CHI-SQUARE-TOL for user \ setting of this parameter; estimated parameter \ uncertainty array (sigmaa) and fitted curve array \ (yfit) are now external inputs; fitting function \ xt is passed to INIT-CURVEFIT; removed MAX_POINTS \ and check on NPTS since relevant arrays are all external. \ \ NOTES specific to this version of curvefit: \ \ (1) INIT-CURVEFIT must be called first, prior to calling CURFIT. \ \ (2) SET-CHI-SQUARE-TOL may be used to set a non-default \ tolerance for the termination of CURFIT: when chi-square \ increases by an amount less than the fraction tolerance, \ CURFIT will not increase lambda and try again. The \ selection of initial parameter increments is important! CR .( CURVEFIT V2.0 14 November 2023 KM ) BEGIN-MODULE BASE @ DECIMAL \ ------ Utilities for working with FSL arrays and matrices : }}size@ ( a -- m | return number of columns in a matrix ) 2 CELLS - @ ; : }}size! ( m a -- | store the number of columns in the matrix, effectively resizing it ) 2 CELLS - ! ; \ --------- End of array and matrix utilities Public: 32 constant MAX_PARAMETERS MAX_PARAMETERS MAX_PARAMETERS FLOAT matrix alpha{{ \ curvature matrix MAX_PARAMETERS MAX_PARAMETERS FLOAT matrix alph2{{ \ modified curv. matrix and inverse MAX_PARAMETERS FLOAT array beta{ \ derivs of chi0^2 w.r.t. params MAX_PARAMETERS FLOAT array bet2{ MAX_PARAMETERS FLOAT array deriv{ \ derivs of f(x) w.r.t. params Private: DEFER functn 0 ptr xa{ \ address of x array (npts) 0 ptr ya{ \ address of y array (npts) 0 ptr yfit{ \ address of yfit array (npts) 0 ptr aa{ \ address of parameter array (nterms) 0 ptr adel{ \ address of parameter increment array (nterms) 0 ptr asig{ \ address of parameter uncertainty array (nterms) \ Default fractional tolerance on chi square for fit termination 1E-5 fconstant DEF_TOL_CHI_SQUARE fvariable tol_chi_square Public: variable npts variable nterms variable nfree fvariable lambda fvariable chisq1 fvariable chisqr \ Arguments to INIT-CURVEFIT are the following: \ \ x array for x values input \ y array for y values input \ yfit array for receiving computed fit curve output \ a parameter array input/output \ deltaa parameter increment array input \ sigmaa estimated parameter uncertainty array output \ nterms number of parameters input \ npts number of points: (x,y) pairs input \ xtfunc execution token for fitting function input \ \ The calling code must set up all of the above arrays and \ define the fitting function prior to calling INIT-CURVEFIT : init-curvefit ( x y yfit a deltaa sigmaa nterms npts xtfunc -- ) IS functn npts ! nterms ! TO asig{ TO adel{ TO aa{ TO yfit{ TO ya{ TO xa{ nterms @ MAX_PARAMETERS > ABORT" Too many parameters for CURFIT" npts @ nterms @ - nfree ! \ Set up sizes of alpha{{ and alph2{{ matrices nterms @ alpha{{ }}size! nterms @ alph2{{ }}size! 0.001e lambda f! DEF_TOL_CHI_SQUARE tol_chi_square f! ; \ Provide user setting of non-default fractional \ chi-square tolerance for termination of CURFIT. : set-chi-square-tol ( F: r -- ) fdup F0> IF tol_chi_square f! ELSE fdrop THEN ; \ Evaluate reduced chi-square for fit to data : fchisq ( -- chi-square) 0.0e0 nfree @ 0> IF npts @ 0 DO \ Accumulate chi-square ya{ I } F@ yfit{ I } F@ F- fsquare F+ LOOP nfree @ s>f F/ \ divide by nfree THEN ; \ Non-analytical derivative routine fvariable aj fvariable delta : fderiv ( n -- | evaluate derivative of function at x_n ) xa{ SWAP } F@ nterms @ 0 DO aa{ I } F@ aj F! adel{ I } F@ FDUP F0= IF FDROP 0.0e0 \ derivative taken to be zero if parameter inc is zero ELSE delta F! aj F@ delta F@ F+ aa{ I } F! FDUP aa{ functn aj F@ delta F@ F- aa{ I } F! FOVER aa{ functn F- delta F@ F/ 2.0e0 F/ aj F@ aa{ I } F! THEN deriv{ I } F! LOOP FDROP ; : curfit ( F: -- chisqr ) \ See Note (2) \ Initialize alpha matrix and beta array to zero. nterms @ 0 DO 0.0e0 beta{ I } F! I 1+ 0 DO 0.0e0 alpha{{ J I }} F! LOOP LOOP \ Compute new derivs and curvature matrix elements. npts @ 0 DO I fderiv nterms @ 0 DO ya{ J } F@ xa{ J } F@ aa{ functn F- deriv{ I } F@ F* beta{ I } F+! I 1+ 0 DO deriv{ J } F@ deriv{ I } F@ F* alpha{{ J I }} F+! LOOP LOOP LOOP nterms @ 0 DO I 1+ 0 DO alpha{{ J I }} F@ alpha{{ I J }} F! LOOP LOOP \ Evaluate chi square at starting point npts @ 0 DO xa{ I } F@ aa{ functn yfit{ I } F! LOOP fchisq chisq1 F! \ Compute modified curvature matrix and invert it \ to find new parameters BEGIN nterms @ 0 DO nterms @ 0 DO alpha{{ J I }} F@ alpha{{ J J }} F@ alpha{{ I I }} F@ F* FDUP F0= IF F2DROP 0.0e0 ELSE FSQRT F/ THEN alph2{{ J I }} F! LOOP lambda F@ 1.0e0 F+ alph2{{ I I }} F! LOOP alph2{{ nterms @ mat^-1 ABORT" Singular Matrix!" nterms @ 0 DO aa{ I } F@ bet2{ I } F! nterms @ 0 DO beta{ I } F@ alph2{{ J I }} F@ F* alpha{{ J J }} F@ alpha{{ I I }} F@ F* FDUP F0= IF F2DROP 0.0e0 ELSE FSQRT F/ THEN bet2{ J } F+! LOOP LOOP \ If chi-square increased, increase lambda and try again npts @ 0 DO xa{ I } F@ bet2{ functn yfit{ I } F! LOOP fchisq chisqr F! chisq1 F@ chisqr F@ f2dup F< >r F- FABS chisq1 F@ F/ tol_chi_square F@ F> r> and WHILE lambda F@ 10.0e0 F* lambda F! REPEAT \ Evaluate parameters and uncertainties nterms @ 0 DO bet2{ I } F@ aa{ I } F! adel{ I } F@ F0= IF 0.0e0 ELSE alph2{{ I I }} F@ alpha{{ I I }} F@ F/ FSQRT THEN asig{ I } F! LOOP lambda F@ 10.0e0 F/ lambda F! chisqr F@ \ return chi-square on stack ; BASE ! END-MODULE