\ Singular Value Decomposition; dealing with (almost) singular sets of \ equations or matrices. \ Forth Scientific Library Algorithm #39 \ Complement to GAUSSJ.xxx, LU.xxx and MARQUARD.xxx \ svdcmp ( 'a m n 'w 'v -- ) \ Given a matrix a[0..m-1,0..n-1] this routine computes its singular value \ decomposition, a = u * w * v^T. The matrix u replaces a . The diagonal \ matrix of singular values w (all positive) is stored in the vector w[0..n-1]. \ The matrix v (not the transpose v^T) is stored as v[0..n-1,0..n-1]. \ svbksb ( 'u 'w 'v m n 'b 'x -- ) \ Solves AX = B for a vector X, where A is specified by the arrays \ u[0..m-1,0..n-1], w[0..n-1], v[0..n-1,0..n-1] as returned by svdcmp. m and \ n are the dimensions of A, and will be equal for square matrices. b[0..m-1] \ is the input right-hand side. x[0..n-1] is the solution vector. No \ input quantities are destroyed, so the routine may be called sequentially \ with different b's. \ svdfit ( 'x 'y 'sig ndata 'a ma 'u 'v 'w -- ) ( F: -- chisq ) \ Given a set of points x[0..ndata-1], y[0..ndata-1] with individual \ standard deviations given by sig[0..ndata-1], use chi-squared minimization \ to determine the coefficients a[0..ma-1] of the fitting function \ y = SUM a[i]*afunc[i](x). Here we solve the fitting equations using singular \ value decomposition of the ndata by ma matrix. Initially, arrays \ u[0..ndata-1,0..ma-1], v[0..ma-1,0..ma-1] and w[0..ma-1] provide workspace; \ when svdfit finishes they define the svd and can be used to obtain the \ covariance matrix. The routine puts computed values for the ma fit \ parameters in a, and returns chi-squared (chisq) on the stack. \ The user supplies a word svdfuncs ( 'afunc ma -- ) ( F: x -- ) that \ returns the ma basic functions evaluated at x in the array afunc[0..ma-1]. \ svdfuncs ( 'afunc ma -- ) ( F: x -- ) \ see sdvfit . The user supplies this routine. It returns the ma basic \ functions evaluated at x in the array afunc[0..ma-1]. \ svdtol ( -- addr ) \ A floating-point variable that controls the cut-off limit of the sdv \ algorithm. Singular values smaller than the product of svdtol and the \ maximum w (see svdcmp) are treated as zero within svdfit . \ A typical value would be 1e-5. \ svdvar ( 'v ma 'w 'cvm -- ) \ To evaluate the covariance matrix cvm[0..ma-1,0..ma-1] of the fit for ma \ parameters obtained by svdfit, call this routine with matrices \ v[0..ma-1,0..ma-1], w[0..ma-1] as set up by svdfit . The user must \ pre-allocate the cvm array. \ This is an ANS Forth program requiring: \ 1. The Floating-Point word sets \ 2. The LOCALS EXT word set (#LOCALS must be 9 or greater) \ 3. [IF] and [THEN] from TOOLS EXT \ 4. Words from fsl_util.xxx \ 5. Memory allocation words from dynmem.xxx \ 6. Uses \ : F0= ( -- bool ) ( F: r -- ) 0e 0e F~ ; \ : F0<> ( -- bool ) ( F: r -- ) F0= 0= ; \ : F> ( -- bool ) ( F: r1 r2 -- ) FSWAP F< ; \ : F2DUP ( F: r1 r2 -- r1 r2 r1 r2 ) FOVER FOVER ; \ : 1/F ( F: r -- 1/r ) 1e FSWAP F/ ; \ : F+! ( addr -- ) ( F: r -- ) DUP F@ F+ F! ; \ : FSQR ( F: r -- r^2 ) FDUP F* ; \ : F1+ ( F: r -- r+1 ) 1e F+ ; \ A random generator for the test words, for instance: \ HERE VALUE seed ( Note: '$' means 'in HEX' ) \ : RANDOM ( -- u ) seed $107465 * $234567 + DUP TO seed ; \ : CHOOSE ( n -- u ) ( 0 <= u < n ) RANDOM UM* NIP ; \ 7. Uses FVALUE ( analogous to VALUE ) \ 8. The test code uses mat* fillmat and transpose from gaussj.xxx \ Note1: the code uses 4 fp stack cells (iForth vsn 1.05) when executing \ the test words. \ Note2: svdfit uses 9 locals, while only 8 might be available in some systems. \ See: 'Numerical recipes in Pascal, The Art of Scientific Computing', \ William H. Press, Brian P. Flannery, Saul A. Teukolsky and William \ T. Vetterling, Chapter 2 (2.9): Singular Value Decomposition, Chapter \ 14 (14.3): General Least Squares, 1989; Cambridge University Press, \ Cambridge, ISBN 0-521-37516-9 \ (c) Copyright 1995 Marcel Hendrix. Permission is granted by the \ author to use this software for any application provided this \ copyright notice is preserved. CR .( SVD V1.1 24 June 1995 MH ) Private: FLOAT DARRAY tmp{ : FSIGN ( F: a b -- c ) \ c = sign(b)*|a| F0< FABS IF FNEGATE THEN ; : pythag ( F: a b -- c ) \ compute sqrt(a^2+b^2) without overflow FABS FSWAP FABS FRAME| a b | a b F> IF b a F/ FSQR F1+ FSQRT a F* |FRAME EXIT THEN b F0= IF 0e ELSE a b F/ FSQR F1+ FSQRT b F* THEN |FRAME ; 0e FVALUE xx 0e FVALUE yy 0e FVALUE zz 0e FVALUE cc 0e FVALUE f \ These could be done with FRAME| 0e FVALUE g 0e FVALUE h FVARIABLE ss FVARIABLE scale FVARIABLE anorm Public: \ Given a matrix a[0..m-1,0..n-1] this routine computes its singular value \ decomposition, a = u * w * v^T. The matrix u replaces a . The diagonal \ matrix of singular values w (all positive) is stored in the vector w[0..n-1]. \ The matrix v (not the transpose v^T) is stored as v[0..n-1,0..n-1]. : svdcmp ( 'a m n 'w 'v -- ) 0 0 0 LOCALS| nm ii ll v{{ w{ n m a{{ | 0e scale F! 0e anorm F! \ Householder reduction to bidiagonal form & tmp{ n }malloc malloc-fail? ABORT" svdcmp :: out of memory" n 0 DO I TO ii I 1+ TO ll scale F@ g F* tmp{ I } F! 0e TO g 0e ss F! 0e scale F! I m < IF m I ?DO a{{ I J }} F@ FABS scale F+! LOOP scale F@ F0<> IF m I ?DO a{{ I J }} DUP F@ scale F@ F/ FDUP F! FSQR ss F+! LOOP a{{ I I }} F@ TO f ss F@ FSQRT f FSIGN FNEGATE TO g f g F* ss F@ F- TO h f g F- a{{ I I }} F! n ll ?DO 0e m ii ?DO a{{ I ii }} F@ a{{ I J }} F@ F* F+ LOOP h F/ ( F: f) m ii ?DO FDUP a{{ I ii }} F@ F* a{{ I J }} F+! LOOP FDROP LOOP m I ?DO a{{ I J }} DUP F@ scale F@ F* F! LOOP THEN THEN scale F@ g F* w{ I } F! 0e TO g 0e ss F! 0e scale F! I m < I n 1- <> AND IF n ll ?DO a{{ J I }} F@ FABS scale F+! LOOP scale F@ F0<> IF n ll ?DO a{{ ii I }} DUP F@ scale F@ F/ FDUP F! FSQR ss F+! LOOP a{{ I ll }} F@ TO f ss F@ FSQRT f FSIGN FNEGATE TO g f g F* ss F@ F- TO h f g F- a{{ I ll }} F! n ll ?DO a{{ ii I }} F@ h F/ tmp{ I } F! LOOP m ll ?DO 0e n ll ?DO a{{ J I }} F@ a{{ ii I }} F@ F* F+ LOOP ( F: ss) n ll ?DO FDUP tmp{ I } F@ F* a{{ J I }} F+! LOOP FDROP LOOP n ll ?DO a{{ J I }} DUP F@ scale F@ F* F! LOOP THEN THEN w{ I } F@ FABS tmp{ I } F@ FABS F+ anorm F@ FMAX anorm F! LOOP \ accumulation of right-hand transforms 0 n 1- DO I TO ii I n 1- < IF g F0<> IF n ll ?DO a{{ J I }} F@ a{{ J ll }} F@ F/ g F/ v{{ I J }} F! LOOP n ll ?DO 0e n ll ?DO a{{ ii I }} F@ v{{ I J }} F@ F* F+ LOOP n ll ?DO FDUP v{{ I ii }} F@ F* v{{ I J }} F+! LOOP FDROP LOOP THEN n ll ?DO 0e v{{ J I }} F! 0e v{{ I J }} F! LOOP THEN 1e v{{ I I }} F! tmp{ I } F@ TO g I TO ll -1 +LOOP \ accumulation of left-hand transforms 0 m n MIN 1- DO I TO ii I 1+ TO ll w{ I } F@ TO g n ll ?DO 0e a{{ J I }} F! LOOP g F0<> IF g 1/F TO g n ll ?DO 0e m ll ?DO a{{ I ii }} F@ a{{ I J }} F@ F* F+ LOOP a{{ J J }} F@ F/ g F* m ii ?DO FDUP a{{ I ii }} F@ F* a{{ I J }} F+! LOOP FDROP LOOP m I ?DO a{{ I J }} DUP F@ g F* F! LOOP ELSE m I ?DO 0e a{{ I J }} F! LOOP THEN 1e a{{ I DUP }} F+! -1 +LOOP \ diagonalization of the bidiagonal form 0 n 1- DO \ loop over singular values I TO ii 31 1 DO \ try 30 times 0 TO ll 1 0 J DO \ test for splitting I 1- TO nm tmp{ I } F@ FABS anorm F@ F+ anorm F@ F= IF DROP I TO ll 0 LEAVE THEN nm -1 > IF w{ nm } F@ FABS anorm F@ F+ anorm F@ F= IF DROP I TO ll 1 LEAVE THEN THEN -1 +LOOP IF \ 1: 0e TO cc 1e ss F! J 1+ ll ?DO tmp{ I } DUP F@ FDUP ss F@ F* TO f cc F* F! f FABS anorm F@ F+ anorm F@ F= IF LEAVE THEN w{ I } F@ TO g f g pythag FDUP TO h w{ I } F! h 1/F TO h g h F* TO cc f h F* FNEGATE ss F! m 0 DO a{{ I nm }} F@ ( y) a{{ I J }} F@ ( y z) FOVER cc F* FOVER ss F@ F* F+ a{{ I nm }} F! cc F* FSWAP ss F@ F* FNEGATE F+ a{{ I J }} F! LOOP LOOP THEN \ 2: w{ J } F@ TO zz ll J = IF zz F0< IF zz FNEGATE w{ J } F! n 0 DO v{{ I ii }} DUP F@ FNEGATE F! LOOP THEN LEAVE THEN I 30 = IF CR ." svdcmp :: doesn't converge" UNLOOP LEAVE THEN w{ ll } F@ TO xx J 1- TO nm w{ nm } F@ TO yy tmp{ nm } F@ TO g tmp{ J } F@ TO h yy zz F- yy zz F+ F* g h F- g h F+ F* F+ h yy F* F2* F/ TO f f 1e pythag TO g g f FSIGN f F+ yy FSWAP F/ h F- h F* xx zz F- xx zz F+ F* F+ xx F/ TO f \ next QR transformation 1e TO cc 1e ss F! nm 1+ ll ?DO I 1+ TO ii tmp{ ii } F@ TO g w{ ii } F@ TO yy ss F@ g F* TO h cc g F* TO g f h pythag FDUP TO zz tmp{ I } F! f zz F/ TO cc h zz F/ ss F! xx cc F* g ss F@ F* F+ TO f xx ss F@ F* FNEGATE g cc F* F+ TO g yy ss F@ F* TO h yy cc F* TO yy n 0 DO v{{ I J }} DUP F@ TO xx v{{ I ii }} DUP F@ TO zz xx ss F@ F* FNEGATE zz cc F* F+ F! xx cc F* zz ss F@ F* F+ F! LOOP f h pythag FDUP TO zz w{ I } F! zz F0<> IF zz 1/F FDUP TO zz FDUP f F* TO cc h F* ss F! THEN cc g F* ss F@ yy F* F+ TO f ss F@ g F* FNEGATE cc yy F* F+ TO xx m 0 DO a{{ I J }} DUP F@ TO yy a{{ I ii }} DUP F@ TO zz yy ss F@ F* FNEGATE zz cc F* F+ F! yy cc F* zz ss F@ F* F+ F! LOOP LOOP 0e tmp{ ll } F! f tmp{ J } F! xx w{ J } F! LOOP -1 +LOOP & tmp{ }free ; \ Solves AX = B for a vector X, where A is specified by the arrays \ u[0..m-1,0..n-1], w[0..n-1], v[0..n-1,0..n-1] as returned by svdcmp. m and \ n are the dimensions of A, and will be equal for square matrices. b[0..m-1] \ is the input right-hand side. x[0..n-1] is the solution vector. No \ input quantities are destroyed, so the routine may be called sequentially \ with different b's. : svbksb ( 'u 'w 'v m n 'b 'x -- ) LOCALS| x{ b{ n m v{{ w{ u{{ | & tmp{ n }malloc malloc-fail? ABORT" svbksb :: out of memory" n 0 DO \ calculate U^T B .. 0e w{ I } F@ F0<> IF m 0 DO u{{ I J }} F@ b{ I } F@ F* F+ LOOP w{ I } F@ F/ THEN tmp{ I } F! LOOP n 0 DO \ .. multiply with V -> answer 0e n 0 DO v{{ J I }} F@ tmp{ I } F@ F* F+ LOOP x{ I } F! LOOP & tmp{ }free ; \ A floating-point variable that controls the cut-off limit of the sdv \ algorithm. Singular values smaller than the product of svdtol and the \ maximum w (see svdcmp) are treated as zero within svdfit . \ A typical value would be 1e-5. FVARIABLE svdtol 1e-32 svdtol F! \ See sdvfit . The user supplies this routine. It returns the ma basic \ functions evaluated at x in the array afunc[0..ma-1]. v: svdfuncs ( 'afunc ma -- ) ( F: x -- ) Private: FLOAT DARRAY afunc{ FLOAT DARRAY b{ Public: \ Given a set of points x[0..ndata-1], y[0..ndata-1] with individual \ standard deviations given by sig[0..ndata-1], use chi-squared minimization \ to determine the coefficients a[0..ma-1] of the fitting function \ y = SUM a[i]*afunc[i](x). Here we solve the fitting equations using singular \ value decomposition of the ndata by ma matrix. Initially, arrays \ u[0..ndata-1,0..ma-1], v[0..ma-1,0..ma-1] and w[0..ma-1] provide workspace; \ when svdfit finishes they define the svd and can be used to obtain the \ covariance matrix. The routine puts computed values for the ma fit \ parameters in a, and returns chi-squared (chisq) on the stack. \ The user supplies a word svdfuncs ( 'afunc ma -- ) ( F: x -- ) that \ returns the ma basic functions evaluated at x in the array afunc[0..ma-1]. : svdfit ( 'x 'y 'sig ndata 'a ma 'u 'v 'w -- ) ( F: -- chisq ) LOCALS| w{ v{{ u{{ ma a{ ndata sig{ y{ x{ | & b{ ndata }malloc malloc-fail? & afunc{ ma }malloc malloc-fail? OR ABORT" svdfit :: out of memory" ndata 0 DO x{ I } F@ afunc{ ma svdfuncs sig{ I } F@ 1/F ma 0 DO FDUP afunc{ I } F@ F* u{{ J I }} F! LOOP y{ I } F@ F* b{ I } F! LOOP u{{ ndata ma w{ v{{ svdcmp 0e ma 0 DO w{ I } F@ FMAX LOOP svdtol F@ F* ( treshold) ma 0 DO w{ I } F@ FOVER F< IF 0e w{ I } F! THEN LOOP FDROP u{{ w{ v{{ ndata ma b{ a{ svbksb 0e ( F: chisq) ndata 0 DO x{ I } F@ afunc{ ma svdfuncs 0e ( F: chisq sum) ma 0 DO a{ I } F@ afunc{ I } F@ F* F+ LOOP y{ I } F@ F- sig{ I } F@ F/ FSQR F+ LOOP ( F: chisq) & afunc{ }free & b{ }free ; \ To evaluate the covariance matrix cvm[0..ma-1,0..ma-1] of the fit for ma \ parameters obtained by svdfit, call this routine with matrices \ v[0..ma-1,0..ma-1], w[0..ma-1] as set up by svdfit . The user must \ pre-allocate the cvm array. : svdvar ( 'v ma 'w 'cvm -- ) 0 LOCALS| ii cvm{{ w{ ma v{{ | & tmp{ ma }malloc malloc-fail? ABORT" svdvar :: out of memory" ma 0 DO w{ I } F@ FDUP F0<> IF FSQR 1/F THEN tmp{ I } F! LOOP ma 0 DO I TO ii I 1+ 0 DO 0e ma 0 DO v{{ ii I }} F@ v{{ J I }} F@ F* tmp{ I } F@ F* F+ LOOP FDUP cvm{{ J I }} F! cvm{{ I J }} F! LOOP LOOP & tmp{ }free ; Reset_Search_Order TEST-CODE? [IF] \ =========================================================== FALSE TO TEST-CODE? \The test code needs mat* fillmat and transpose from gaussj.frt S" gaussj.frt" INCLUDED TRUE TO TEST-CODE? \ define suitable svdfuncs \ Simple polynomial of degree NP-1 with coefficients in the array p{ : fpoly ( 'a n -- ) ( F: x -- ) LOCALS| np p{ | 1e p{ 0 } F! ( F: x) np 1 DO FDUP p{ I 1- } F@ F* p{ I } F! LOOP FDROP ; \ Simple polynomial of degree 2NP-1 with coefficients in the array p{ : fpoly[odd] ( 'a n -- ) ( F: x -- ) LOCALS| np p{ | FDUP p{ 0 } F! FSQR ( F: x^2) np 1 DO FDUP p{ I 1- } F@ F* p{ I } F! LOOP FDROP ; \ Expansion with nl Legendre polynomials pl : fleg ( 'a n -- ) ( F: x -- ) LOCALS| nl pl{ | FRAME| a | 1e pl{ 0 } F! a pl{ 1 } F! nl 3 < IF |FRAME EXIT THEN a F2* &b F! a &c F! 1e &d F! nl 2 DO d &e F! b &c F+! 1e &d F+! pl{ I 1- } F@ c F* pl{ I 2- } F@ e F* F- d F/ pl{ I } F! LOOP |FRAME ; & fleg defines svdfuncs 10 VALUE #dp \ don't make it bigger than 60 (allocation below!) 3 VALUE #vr 60 FLOAT ARRAY x{ 60 FLOAT ARRAY y{ 60 FLOAT ARRAY sig{ 60 FLOAT ARRAY a{ FLOAT DARRAY w{ FLOAT DMATRIX u{{ FLOAT DMATRIX v{{ FLOAT DMATRIX a{{ FLOAT DMATRIX w{{ FLOAT DMATRIX svdcovar{{ FLOAT DMATRIX mult{{ FLOAT DMATRIX vT{{ FLOAT DMATRIX a'{{ : alloc #dp 60 > #vr 20 > OR ABORT" svd :: problem too large for x{ y{ sig{ a{" & w{ #vr }malloc malloc-fail? & u{{ #dp #vr }}malloc malloc-fail? OR & v{{ #vr #vr }}malloc malloc-fail? OR & a{{ #dp #vr }}malloc malloc-fail? OR & w{{ #vr #vr }}malloc malloc-fail? OR & svdcovar{{ #vr #vr }}malloc malloc-fail? OR ABORT" svd :: out of memory" ; : dealloc & svdcovar{{ }}free & w{{ }}free & a{{ }}free & v{{ }}free & u{{ }}free & w{ }free ; \ Assemble a random matrix a{{ , decompose it into u{{ v{{ and w{{ , then \ reassemble it (using matrix multiplication) to check if the same matrix \ emerges again. (Normally one inspects the w{ matrix for vanishingly small \ elements, zeroing them before attempting the reassembly). \ Because v{{ and u{{ are orthonormal we also check if v{{ * vT{{ = I, \ u{{ * uT{{ = I (unit matrix). : (TEST-SVD) a{{ u{{ #dp #vr }}fcopy u{{ #dp #vr w{ v{{ svdcmp w{{ #vr #vr 0e fillmat #vr 0 DO w{ I } F@ w{{ I I }} F! LOOP u{{ #dp #vr w{{ #vr #vr & mult{{ mat* v{{ #vr #vr & vT{{ transpose mult{{ #dp #vr vT{{ #vr #vr & a'{{ mat* CR ." # datapoints = " #dp . ." # variables = " #vr . CR ." U:" #dp #vr u{{ CR }}fprint CR ." W:" #vr #vr w{{ CR }}fprint CR ." V^T:" #vr #vr vT{{ CR }}fprint CR ." A should be equal to A'" CR ." A:" #dp #vr a{{ CR }}fprint CR ." A':" #dp #vr a'{{ CR }}fprint & a'{{ }}free & mult{{ }}free vT{{ #vr #vr v{{ #vr #vr & mult{{ mat* CR ." V^T * V (= I) : " #vr #vr mult{{ CR }}fprint & mult{{ }}free & vT{{ }}free u{{ #dp #vr & vT{{ transpose vT{{ #vr #dp u{{ #dp #vr & mult{{ mat* CR ." U^T * U (= I) : " #vr #vr mult{{ CR }}fprint & mult{{ }}free & vT{{ }}free ; : TEST-SVDCMP ( -- ) 60 TO #dp 4 TO #vr alloc #dp 0 DO #vr 0 DO 10 CHOOSE S>F a{{ J I }} F! LOOP LOOP (TEST-SVD) dealloc ; \ Deliberately make the array singular by zeroing a column. : TEST-SINGULAR ( -- ) 4 TO #dp 4 TO #vr alloc #dp 0 DO #vr 0 DO 10 CHOOSE S>F a{{ J I }} F! LOOP LOOP #dp 0 DO 0e a{{ I 1 }} F! LOOP (TEST-SVD) dealloc ; : PROBLEM1 ( -- ) alloc & fpoly defines svdfuncs #dp 0 DO I S>F 5e #dp S>F F/ F* FDUP x{ I } F! FDUP 2.33e F* FSWAP FSQR 1.34567e F* F+ y{ I } F! 1e-2 sig{ I } F! LOOP x{ y{ sig{ #dp a{ #vr u{{ v{{ w{ svdfit CR ." # datapoints = " #dp . ." # variables = " #vr . CR ." Fitter returns chi-squared = " F. CR ." Coefficients (polynomial, expected 0, 2.33, 1.34567): " #vr a{ }fprint v{{ #vr w{ svdcovar{{ svdvar CR CR ." covariances: " #vr #vr svdcovar{{ CR }}fprint dealloc ; CR .( Try: PROBLEM1 ) \ Copied from polys.xxx : Le_n ( n -- ) ( F: x -- z ) \ nth order Legendre Polynomial FDUP 1e FRAME| a b c | DUP 0= IF DROP 1e &b F! ELSE DUP 1 > IF 1 DO c b F* FDUP a F- I S>F I 1+ S>F F/ F* F+ b &a F! &b F! LOOP ELSE DROP THEN THEN b |FRAME ; : }Horner(odd) ( &a n -- ) ( F: x -- y[x] ) FDUP FSQR FRAME| a | 0e 0 SWAP 1- DO a F* DUP I } F@ F+ -1 +LOOP DROP F* |FRAME ; : PROBLEM2 ( -- ) 0e FRAME| a | alloc & fleg defines svdfuncs #dp 0 DO I S>F 5e #dp S>F F/ F* FDUP x{ I } F! FDUP 0 Le_n 3.4356e F* &a F! FDUP 1 Le_n -3.4356e F* &a F+! 2 Le_n 8e F* &a F+! a y{ I } F! 1e-2 sig{ I } F! LOOP x{ y{ sig{ #dp a{ #vr u{{ v{{ w{ svdfit CR ." # datapoints = " #dp . ." # variables = " #vr . CR ." Fitter returns chi-squared = " F. CR ." Coefficients (Legendre, expected 3.4356, -3.4356, 8): " #vr a{ }fprint v{{ #vr w{ svdcovar{{ svdvar CR CR ." covariances: " #vr #vr svdcovar{{ CR }}fprint |FRAME dealloc ; CR .( Try: PROBLEM2 ) : )FIT-SINE ( xt -- ) #vr LOCALS| ma | alloc defines svdfuncs #dp 0 DO I S>F PI/2 #dp S>F F/ F* FDUP x{ I } F! FSIN y{ I } F! 1e-2 sig{ I } F! \ don't make this too small! LOOP x{ y{ sig{ #dp a{ ma u{{ v{{ w{ svdfit CR ." # datapoints = " #dp . ." # variables = " #vr . CR ." Fitter returns chi-squared = " F. CR ." Coefficients: " ma a{ }fprint CR ." w array: " ma w{ }fprint dealloc ; FALSE [IF] ( For those without XY-PLOT) : XY-PLOT ( 'x 'y n -- ) LOCALS| n y{ x{ | n 0 DO CR I 4 .R x{ I } F@ F. y{ I } F@ F. LOOP ; [THEN] : TEST-SINE #dp 0 DO x{ I } F@ a{ #vr }Horner(odd) y{ I } F@ F- sig{ I } F! LOOP x{ sig{ #dp XY-PLOT ; CR .( Try: 1e-32 svdtol F! use( fpoly[odd] ) CHAR ) EMIT .( FIT-SINE ) 0 [IF] ( Example output ) FORTH> PROBLEM1 # datapoints = 4 # variables = 4 Fitter returns chi-squared = 9.786174E-32 Coefficients (polynomial, expected 0, 2.33, 1.34567): -1.504331E-0018 2.330000E+0000 1.345670E+0000 -4.201283E-0019 covariances: 1.000000E-0004 -1.466667E-0004 6.400000E-0005 -8.533333E-0006 -1.466667E-0004 9.422222E-0004 -6.400000E-0004 1.069511E-0004 6.400000E-0005 -6.400000E-0004 4.710400E-0004 -8.192000E-0005 -8.533333E-0006 1.069511E-0004 -8.192000E-0005 1.456356E-0005 ok FORTH> PROBLEM2 # datapoints = 4 # variables = 4 Fitter returns chi-squared = 1.833442E-29 Coefficients (Legendre, expected 3.4356, -3.4356, 8): 3.435600E+0000 -3.435600E+0000 8.000000E+0000 2.222614E-0018 covariances: 1.950044E-0004 -3.815040E-0004 1.473422E-0004 -1.433600E-0005 -3.815040E-0004 1.075806E-0003 -4.594347E-0004 4.627570E-0005 1.473422E-0004 -4.594347E-0004 2.093511E-0004 -2.184533E-0005 -1.433600E-0005 4.627570E-0005 -2.184533E-0005 2.330169E-0006 ok FORTH> TEST-SINGULAR # datapoints = 4 # variables = 4 U: -6.272218E-0001 3.308627E-0001 2.777152E-0001 6.480717E-0001 -4.289447E-0001 5.798998E-0001 -4.710812E-0002 -6.910162E-0001 -3.505253E-0001 -5.875159E-0001 6.553417E-0001 -3.201318E-0001 -5.474772E-0001 -4.572430E-0001 -7.008440E-0001 3.904046E-0003 W: 1.588664E+0001 0.000000E+0000 0.000000E+0000 0.000000E+0000 0.000000E+0000 7.625206E+0000 0.000000E+0000 0.000000E+0000 0.000000E+0000 0.000000E+0000 2.114470E+0000 0.000000E+0000 0.000000E+0000 0.000000E+0000 0.000000E+0000 2.098447E-0020 V^T: -7.312427E-0001 6.737937E-0021 -4.896171E-0001 -4.749307E-0001 6.820110E-0001 3.122598E-0020 -5.125054E-0001 -5.217271E-0001 -1.204196E-0002 -9.774396E-0020 7.054171E-0001 -7.086901E-0001 1.754642E-0020 -1.000000E+0000 6.331127E-0021 1.103536E-0019 A should be equal to A' A: 9.000000E+0000 0.000000E+0000 4.000000E+0000 3.000000E+0000 8.000000E+0000 0.000000E+0000 1.000000E+0000 1.000000E+0000 1.000000E+0000 0.000000E+0000 6.000000E+0000 4.000000E+0000 4.000000E+0000 0.000000E+0000 5.000000E+0000 7.000000E+0000 A': 9.000000E+0000 -5.935658E-0020 4.000000E+0000 3.000000E+0000 8.000000E+0000 1.163979E-0019 1.000000E+0000 1.000000E+0000 1.000000E+0000 -3.061375E-0019 6.000000E+0000 4.000000E+0000 4.000000E+0000 -2.270917E-0020 5.000000E+0000 7.000000E+0000 V^T * V (= I) : 1.000000E+0000 0.000000E+0000 0.000000E+0000 -7.507878E-0020 0.000000E+0000 1.000000E+0000 -8.131516E-0020 -8.007833E-0020 0.000000E+0000 -8.131516E-0020 1.000000E+0000 2.379224E-0020 -7.507878E-0020 -8.007833E-0020 2.379224E-0020 1.000000E+0000 U^T * U (= I) : 1.000000E+0000 1.355253E-0019 0.000000E+0000 2.202286E-0020 1.355253E-0019 1.000000E+0000 -2.710505E-0020 -7.729176E-0021 0.000000E+0000 -2.710505E-0020 1.000000E+0000 1.103260E-0019 2.202286E-0020 -7.729176E-0021 1.103260E-0019 1.000000E+0000 ok FORTH> TEST-SVDCMP # datapoints = 60 # variables = 4 U: -1.569359E-0001 1.048297E-0001 1.651167E-0001 -1.356633E-0001 -8.382130E-0002 -1.074453E-0001 5.050485E-0002 -2.548267E-0001 -1.153018E-0001 -8.700374E-0003 -6.906958E-0002 2.180546E-0001 -1.194084E-0001 -7.708235E-0003 1.779695E-0002 -2.687051E-0001 -1.637876E-0001 7.566673E-0002 -9.806892E-0002 1.641742E-0001 -9.037077E-0002 -7.742788E-0002 -7.963653E-0002 -3.007598E-0001 -1.274799E-0001 2.062557E-0002 -2.298035E-0001 1.701644E-0002 -1.554766E-0001 1.439338E-0002 1.255737E-0001 -2.898255E-0002 -1.042211E-0001 1.177421E-0001 1.730752E-0001 -2.635276E-0003 -7.693363E-0002 -1.217853E-0001 -2.224370E-0001 -4.490494E-0002 -7.773425E-0002 -1.363193E-0001 -1.635911E-0001 9.538559E-0002 -1.520565E-0001 -1.251574E-0001 1.027759E-0001 -7.919962E-0002 -1.619460E-0001 1.424888E-0002 1.010289E-0001 -5.416997E-0002 -1.819758E-0001 1.773644E-0002 -1.872887E-0001 1.027279E-0001 -1.556368E-0001 7.471726E-0002 -8.561937E-0002 -7.047381E-0002 -1.366880E-0001 1.622625E-0001 -1.080438E-0001 -1.596908E-0001 -9.909171E-0002 2.206263E-0001 5.481048E-0002 -2.268222E-0002 -1.882048E-0001 -7.289388E-0002 1.049561E-0001 1.397773E-0001 -1.043918E-0001 -1.877449E-0001 -6.151704E-0003 1.249489E-0001 -9.771255E-0002 1.603519E-0001 -9.032905E-0002 6.325289E-0002 -8.780260E-0002 -1.536802E-0002 1.081908E-0001 2.824568E-0002 -7.625449E-0002 -2.630692E-0001 -6.361230E-0003 1.920888E-0001 -1.581958E-0001 1.774075E-0001 1.279036E-0001 -1.601520E-0001 -1.706361E-0001 1.324720E-0001 -1.130762E-0001 -5.046841E-0003 -1.489476E-0001 2.068957E-0002 5.894225E-0002 2.692806E-0002 -1.127993E-0001 -3.157399E-0002 5.676203E-0002 -1.920489E-0001 -1.124919E-0001 1.639345E-0001 2.230883E-0003 2.008943E-0001 -1.063421E-0001 2.271805E-0001 -5.486301E-0002 1.127713E-0001 -1.133810E-0001 -4.273389E-0002 -1.051227E-0001 5.078401E-0002 -1.062498E-0001 -6.159142E-0002 1.869034E-0001 -1.461159E-0001 -1.509681E-0001 6.365238E-0002 1.333734E-0001 1.531027E-0001 -1.832953E-0001 8.416251E-0002 -1.333254E-0001 4.751595E-0002 -8.310302E-0002 -2.062639E-0001 -2.136852E-0002 2.286771E-0002 -1.140592E-0001 4.100377E-0002 1.954897E-0002 -2.165376E-0001 -1.659561E-0001 -7.965655E-0002 -9.866098E-0002 -3.077518E-0002 -9.065248E-0002 -1.729219E-0001 -1.590793E-0002 3.503319E-0002 -1.503677E-0001 1.535912E-0001 -1.894639E-0001 -3.905193E-0002 -1.442480E-0001 -1.702064E-0001 -7.061708E-0002 2.155009E-0002 -8.472160E-0002 -1.130500E-0001 1.477288E-0001 -1.556322E-0001 -1.572576E-0001 2.441665E-0002 1.760470E-0001 1.897569E-0001 -1.347974E-0001 -2.588146E-0001 3.184380E-0002 -1.782364E-0001 -1.046517E-0001 -1.184917E-0001 -1.789667E-0001 4.236163E-0002 -1.308591E-0001 8.754915E-0002 1.865400E-0001 4.727840E-0002 -1.007208E-0001 -5.197065E-0002 2.558739E-0001 -3.210672E-0002 -1.244894E-0001 9.662294E-0002 2.494629E-0001 3.136981E-0002 -1.129512E-0001 2.510462E-0001 -9.382810E-0002 3.611503E-0002 -1.256300E-0001 1.815041E-0001 2.989725E-0002 6.832746E-0002 -5.077670E-0002 -1.692278E-0001 -9.541712E-0002 1.587824E-0001 -7.971243E-0002 3.507709E-0002 -1.289308E-0001 -2.067976E-0001 -9.903209E-0002 2.267780E-0001 -3.636582E-0002 8.040970E-0003 -1.613870E-0001 -8.179214E-0002 1.587098E-0001 1.517053E-0001 -7.184342E-0002 -6.136642E-0002 -5.275260E-0002 -1.056526E-0001 -1.248195E-0001 -2.060865E-0001 1.997903E-0001 8.713479E-0002 -1.100811E-0001 -1.370418E-0001 -2.863153E-0001 -3.055150E-0002 -1.595864E-0001 -7.058276E-0002 -3.573810E-0002 -4.668376E-0002 -1.800746E-0001 -3.752973E-0002 -7.936720E-0002 -8.489312E-0002 -1.915335E-0001 -1.645687E-0001 2.877294E-0002 -3.856478E-0002 -1.213897E-0001 7.771987E-0002 -1.957211E-0001 -1.018298E-0001 -1.356297E-0001 -3.597123E-0002 9.849439E-0002 2.213365E-0001 -8.740386E-0002 -5.108517E-0002 -4.532132E-0002 1.926297E-0001 W: 7.448449E+0001 0.000000E+0000 0.000000E+0000 0.000000E+0000 0.000000E+0000 2.278041E+0001 0.000000E+0000 0.000000E+0000 0.000000E+0000 0.000000E+0000 2.554021E+0001 0.000000E+0000 0.000000E+0000 0.000000E+0000 0.000000E+0000 2.430660E+0001 V^T: -4.788835E-0001 -4.818671E-0001 -5.623180E-0001 -4.714585E-0001 -3.468425E-0001 -3.291828E-0003 7.595450E-0001 -5.502553E-0001 7.215984E-0001 -6.268806E-0001 1.394648E-0001 -2.585847E-0001 -3.600926E-0001 -6.122205E-0001 2.957015E-0001 6.388114E-0001 A should be equal to A' A: 9.000000E+0000 5.000000E+0000 8.000000E+0000 1.000000E+0000 7.000000E+0000 6.000000E+0000 0.000000E+0000 0.000000E+0000 1.000000E+0000 2.000000E+0000 6.000000E+0000 8.000000E+0000 7.000000E+0000 8.000000E+0000 3.000000E+0000 0.000000E+0000 2.000000E+0000 5.000000E+0000 9.000000E+0000 8.000000E+0000 5.000000E+0000 9.000000E+0000 0.000000E+0000 0.000000E+0000 0.000000E+0000 8.000000E+0000 5.000000E+0000 6.000000E+0000 8.000000E+0000 4.000000E+0000 7.000000E+0000 4.000000E+0000 6.000000E+0000 1.000000E+0000 7.000000E+0000 1.000000E+0000 0.000000E+0000 7.000000E+0000 0.000000E+0000 5.000000E+0000 0.000000E+0000 4.000000E+0000 1.000000E+0000 7.000000E+0000 9.000000E+0000 5.000000E+0000 4.000000E+0000 5.000000E+0000 8.000000E+0000 5.000000E+0000 7.000000E+0000 4.000000E+0000 2.000000E+0000 8.000000E+0000 8.000000E+0000 9.000000E+0000 4.000000E+0000 8.000000E+0000 7.000000E+0000 4.000000E+0000 3.000000E+0000 9.000000E+0000 7.000000E+0000 1.000000E+0000 3.000000E+0000 3.000000E+0000 8.000000E+0000 0.000000E+0000 8.000000E+0000 3.000000E+0000 8.000000E+0000 9.000000E+0000 4.000000E+0000 2.000000E+0000 2.000000E+0000 8.000000E+0000 0.000000E+0000 4.000000E+0000 7.000000E+0000 3.000000E+0000 5.000000E+0000 1.000000E+0000 4.000000E+0000 3.000000E+0000 3.000000E+0000 0.000000E+0000 0.000000E+0000 9.000000E+0000 8.000000E+0000 6.000000E+0000 9.000000E+0000 0.000000E+0000 3.000000E+0000 8.000000E+0000 9.000000E+0000 5.000000E+0000 6.000000E+0000 4.000000E+0000 7.000000E+0000 5.000000E+0000 7.000000E+0000 6.000000E+0000 3.000000E+0000 1.000000E+0000 1.000000E+0000 1.000000E+0000 9.000000E+0000 5.000000E+0000 0.000000E+0000 3.000000E+0000 9.000000E+0000 3.000000E+0000 2.000000E+0000 5.000000E+0000 4.000000E+0000 6.000000E+0000 9.000000E+0000 3.000000E+0000 3.000000E+0000 1.000000E+0000 6.000000E+0000 1.000000E+0000 9.000000E+0000 6.000000E+0000 3.000000E+0000 8.000000E+0000 9.000000E+0000 7.000000E+0000 4.000000E+0000 3.000000E+0000 0.000000E+0000 6.000000E+0000 6.000000E+0000 7.000000E+0000 4.000000E+0000 0.000000E+0000 5.000000E+0000 8.000000E+0000 5.000000E+0000 7.000000E+0000 4.000000E+0000 3.000000E+0000 1.000000E+0000 6.000000E+0000 1.000000E+0000 9.000000E+0000 8.000000E+0000 4.000000E+0000 5.000000E+0000 6.000000E+0000 3.000000E+0000 8.000000E+0000 8.000000E+0000 3.000000E+0000 1.000000E+0000 1.000000E+0000 7.000000E+0000 0.000000E+0000 9.000000E+0000 7.000000E+0000 9.000000E+0000 7.000000E+0000 0.000000E+0000 5.000000E+0000 1.000000E+0000 6.000000E+0000 2.000000E+0000 7.000000E+0000 7.000000E+0000 1.000000E+0000 8.000000E+0000 3.000000E+0000 9.000000E+0000 0.000000E+0000 4.000000E+0000 2.000000E+0000 8.000000E+0000 0.000000E+0000 8.000000E+0000 2.000000E+0000 0.000000E+0000 5.000000E+0000 9.000000E+0000 2.000000E+0000 3.000000E+0000 3.000000E+0000 9.000000E+0000 3.000000E+0000 0.000000E+0000 1.000000E+0000 0.000000E+0000 7.000000E+0000 2.000000E+0000 8.000000E+0000 2.000000E+0000 0.000000E+0000 1.000000E+0000 4.000000E+0000 8.000000E+0000 1.000000E+0000 8.000000E+0000 1.000000E+0000 7.000000E+0000 8.000000E+0000 3.000000E+0000 5.000000E+0000 1.000000E+0000 2.000000E+0000 9.000000E+0000 0.000000E+0000 3.000000E+0000 7.000000E+0000 0.000000E+0000 9.000000E+0000 1.000000E+0000 7.000000E+0000 6.000000E+0000 7.000000E+0000 5.000000E+0000 6.000000E+0000 6.000000E+0000 9.000000E+0000 6.000000E+0000 6.000000E+0000 9.000000E+0000 7.000000E+0000 5.000000E+0000 8.000000E+0000 1.000000E+0000 9.000000E+0000 5.000000E+0000 3.000000E+0000 5.000000E+0000 0.000000E+0000 7.000000E+0000 8.000000E+0000 1.000000E+0000 1.000000E+0000 4.000000E+0000 7.000000E+0000 A': 9.000000E+0000 5.000000E+0000 8.000000E+0000 1.000000E+0000 7.000000E+0000 6.000000E+0000 5.421011E-0019 -2.818926E-0018 1.000000E+0000 2.000000E+0000 6.000000E+0000 8.000000E+0000 7.000000E+0000 8.000000E+0000 3.000000E+0000 -2.168404E-0018 2.000000E+0000 5.000000E+0000 9.000000E+0000 8.000000E+0000 5.000000E+0000 9.000000E+0000 1.301043E-0018 -2.168404E-0018 1.436568E-0018 8.000000E+0000 5.000000E+0000 6.000000E+0000 8.000000E+0000 4.000000E+0000 7.000000E+0000 4.000000E+0000 6.000000E+0000 1.000000E+0000 7.000000E+0000 1.000000E+0000 1.382358E-0018 7.000000E+0000 2.331035E-0018 5.000000E+0000 5.421011E-0020 4.000000E+0000 1.000000E+0000 7.000000E+0000 9.000000E+0000 5.000000E+0000 4.000000E+0000 5.000000E+0000 8.000000E+0000 5.000000E+0000 7.000000E+0000 4.000000E+0000 2.000000E+0000 8.000000E+0000 8.000000E+0000 9.000000E+0000 4.000000E+0000 8.000000E+0000 7.000000E+0000 4.000000E+0000 3.000000E+0000 9.000000E+0000 7.000000E+0000 1.000000E+0000 3.000000E+0000 3.000000E+0000 8.000000E+0000 1.382358E-0018 8.000000E+0000 3.000000E+0000 8.000000E+0000 9.000000E+0000 4.000000E+0000 2.000000E+0000 2.000000E+0000 8.000000E+0000 -1.680513E-0018 4.000000E+0000 7.000000E+0000 3.000000E+0000 5.000000E+0000 1.000000E+0000 4.000000E+0000 3.000000E+0000 3.000000E+0000 -1.734723E-0018 1.301043E-0018 9.000000E+0000 8.000000E+0000 6.000000E+0000 9.000000E+0000 1.517883E-0018 3.000000E+0000 8.000000E+0000 9.000000E+0000 5.000000E+0000 6.000000E+0000 4.000000E+0000 7.000000E+0000 5.000000E+0000 7.000000E+0000 6.000000E+0000 3.000000E+0000 1.000000E+0000 1.000000E+0000 1.000000E+0000 9.000000E+0000 5.000000E+0000 -1.572093E-0018 3.000000E+0000 9.000000E+0000 3.000000E+0000 2.000000E+0000 5.000000E+0000 4.000000E+0000 6.000000E+0000 9.000000E+0000 3.000000E+0000 3.000000E+0000 1.000000E+0000 6.000000E+0000 1.000000E+0000 9.000000E+0000 6.000000E+0000 3.000000E+0000 8.000000E+0000 9.000000E+0000 7.000000E+0000 4.000000E+0000 3.000000E+0000 -1.897354E-0019 6.000000E+0000 6.000000E+0000 7.000000E+0000 4.000000E+0000 -1.517883E-0018 5.000000E+0000 8.000000E+0000 5.000000E+0000 7.000000E+0000 4.000000E+0000 3.000000E+0000 1.000000E+0000 6.000000E+0000 1.000000E+0000 9.000000E+0000 8.000000E+0000 4.000000E+0000 5.000000E+0000 6.000000E+0000 3.000000E+0000 8.000000E+0000 8.000000E+0000 3.000000E+0000 1.000000E+0000 1.000000E+0000 7.000000E+0000 -8.673617E-0019 9.000000E+0000 7.000000E+0000 9.000000E+0000 7.000000E+0000 1.734723E-0018 5.000000E+0000 1.000000E+0000 6.000000E+0000 2.000000E+0000 7.000000E+0000 7.000000E+0000 1.000000E+0000 8.000000E+0000 3.000000E+0000 9.000000E+0000 4.065758E-0019 4.000000E+0000 2.000000E+0000 8.000000E+0000 2.710505E-0019 8.000000E+0000 2.000000E+0000 -1.951564E-0018 5.000000E+0000 9.000000E+0000 2.000000E+0000 3.000000E+0000 3.000000E+0000 9.000000E+0000 3.000000E+0000 5.421011E-0019 1.000000E+0000 6.505213E-0019 7.000000E+0000 2.000000E+0000 8.000000E+0000 2.000000E+0000 -1.084202E-0018 1.000000E+0000 4.000000E+0000 8.000000E+0000 1.000000E+0000 8.000000E+0000 1.000000E+0000 7.000000E+0000 8.000000E+0000 3.000000E+0000 5.000000E+0000 1.000000E+0000 2.000000E+0000 9.000000E+0000 0.000000E+0000 3.000000E+0000 7.000000E+0000 1.951564E-0018 9.000000E+0000 1.000000E+0000 7.000000E+0000 6.000000E+0000 7.000000E+0000 5.000000E+0000 6.000000E+0000 6.000000E+0000 9.000000E+0000 6.000000E+0000 6.000000E+0000 9.000000E+0000 7.000000E+0000 5.000000E+0000 8.000000E+0000 1.000000E+0000 9.000000E+0000 5.000000E+0000 3.000000E+0000 5.000000E+0000 4.336809E-0019 7.000000E+0000 8.000000E+0000 1.000000E+0000 1.000000E+0000 4.000000E+0000 7.000000E+0000 V^T * V (= I) : 1.000000E+0000 2.710505E-0020 3.388132E-0020 5.421011E-0020 2.710505E-0020 1.000000E+0000 -5.421011E-0020 2.710505E-0020 3.388132E-0020 -5.421011E-0020 1.000000E+0000 8.131516E-0020 5.421011E-0020 2.710505E-0020 8.131516E-0020 1.000000E+0000 U^T * U (= I) : 1.000000E+0000 7.199780E-0021 -7.580945E-0020 7.115077E-0020 7.199780E-0021 1.000000E+0000 -8.533857E-0020 -2.456396E-0020 -7.580945E-0020 -8.533857E-0020 1.000000E+0000 -2.795209E-0020 7.115077E-0020 -2.456396E-0020 -2.795209E-0020 1.000000E+0000 ok FORTH> 1e-32 svdtol F! use( fpoly[odd] )FIT-SINE # datapoints = 60 # variables = 4 Fitter returns chi-squared = 6.713799E-8 Coefficients: 9.999977E-0001 -1.666525E-0001 8.310364E-0003 -1.847190E-0004 w array: 5.041953E+0003 6.199818E+0002 1.661439E+0002 2.616607E+0001 ok FORTH> 8 TO #vr ok FORTH> use( fpoly[odd] )FIT-SINE # datapoints = 60 # variables = 8 Fitter returns chi-squared = 1.227569E-27 Coefficients: 1.000000E+0000 -1.666667E-0001 8.333333E-0003 -1.984127E-0004 2.755732E-0006 -2.505187E-0008 1.604790E-0010 -7.372557E-0013 w array: 1.184889E+0005 4.929634E+0003 7.432629E+0002 2.546104E+0002 6.911026E+0001 1.460875E+0001 2.325068E+0000 2.445736E-0001 ok [THEN] [THEN] ( * End of File * )