\ Nelder and Mead's Downhill Simplex Method. \ Forth Scientific Library Algorithm #40 \ amoeba ( 'p 'y xt -- steps ) ( F: tol -- ) \ Multidimensional minimization of the function func(x) where x[0..ndim-1] \ is a vector in ndim dimensions, by the downhill simplex method of \ Nelder and Mead. Set the value NDIM before calling this function. The matrix \ p[0..ndim][0..ndim-1] is input. Its ndim+1 rows are ndim-dimensional vectors \ which are the vertices of the starting simplex. Also input is the vector \ y[0..ndim], whose components must be pre-initialized to the values of func \ evaluated at the ndim+1 vertices (rows) of p; and ftol the fractional \ convergence tolerance to be achieved in the function value (*). \ On output, p and y will have been reset to ndim+1 new points all within ftol \ of a minimum function value, and nfunc gives the number of function \ evaluations taken (nfunc is made -1 to signal non-convergence). \ (*) NB fractional: When the minimum is 0.0 you have a problem! \ This is an ANS Forth program requiring: \ 1. The Floating-Point word sets \ 2. Uses FSL words from fsl_util.xxx \ 3. Uses : F> ( -- bool ) ( F: r1 r2 -- ) FSWAP F< ; \ : F+! ( addr -- ) ( F: r -- ) DUP F@ F+ F! ; \ : FSQR ( F: r -- r^2 ) FDUP F* ; \ : F>= ( F: r1 r2 -- ) ( -- bool ) F< 0= ; \ : F<= ( F: r1 r2 -- ) ( -- bool ) F> 0= ; \ Note: the code uses 5 fp stack cells (iForth vsn 1.05) when executing \ the test words. \ See: 'Numerical Recipes in Pascal, The Art of Scientific Computing', \ William H. Press, Brian P. Flannery, Saul A. Teukolsky and William \ T. Vetterling, Chapter 10 (10.4) 'Downhill Simplex Method in Multidimensions' \ 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 .( AMOEBA V1.1 8 Dec 1995 MH ) Public: 0 VALUE NDIM \ Dimensions of the problem 5000 VALUE nfuncmax \ Change when convergence is very bad. v: TEST-HOOK ( -- ) \ To check inner workings of AMOEBA Private: v: func ( 'array -- ) ( F: -- r ) FLOAT DARRAY ptrY{ \ array 0..ndim-1 scratch FLOAT DARRAY psum{ \ array 0..ndim-1 scratch FLOAT DMATRIX p{{ \ input parameters FLOAT DARRAY y{ 0 VALUE nfunc \ function evaluations used 0 VALUE inhi 0 VALUE ilo 0 VALUE ihi \ Extrapolate by a factor fac through the face of the simplex across from \ the high point, try it, and replace the high point if the new point is \ better. : amoebatry ( -- ) ( F: fac -- ytry ) 0e 0e 0e FRAME| a b c | \ a=ytry 1e FOVER F- NDIM S>F F/ &b F! \ b=fac1 b FSWAP F- &c F! \ c=fac2 NDIM 0 DO psum{ I } F@ b F* p{{ ihi I }} F@ c F* F- ptrY{ I } F! LOOP ptrY{ func &a F! nfunc 1+ TO nfunc a y{ ihi } F@ F< IF a y{ ihi } F! NDIM 0 DO ptrY{ I } F@ FDUP p{{ ihi I }} DUP F@ F- psum{ I } F+! F! LOOP THEN a |FRAME ; : reflect -1.0e amoebatry ; ( F: -- ytry ) : halve 0.5e amoebatry ; ( F: -- ytry ) : double 2.0e amoebatry ; ( F: -- ytry ) \ Look for the worst, the second best and the best function results. : order'm ( -- ) 0 TO ilo y{ 0 } F@ y{ 1 } F@ F> IF 0 TO ihi 1 TO inhi ELSE 1 TO ihi 0 TO inhi THEN NDIM 1+ 0 DO y{ I } F@ y{ ilo } F@ F< IF I TO ilo THEN y{ I } F@ y{ ihi } F@ F> IF ihi TO inhi I TO ihi ELSE y{ I } F@ y{ inhi } F@ F> IF I ihi <> IF I TO inhi THEN THEN THEN LOOP ; \ Compute the geometric mean of the vertices. : geomean ( -- ) NDIM 0 DO 0e NDIM 1+ 0 DO p{{ I J }} F@ F+ LOOP psum{ I } F! LOOP ; \ Compute the fractional range from highest to lowest : fract-range ( F: -- r ) y{ ihi } F@ y{ ilo } F@ F- FABS F2* y{ ihi } F@ FABS y{ ilo } F@ FABS F+ F/ ; : store.p ( -- ) NDIM 1+ 0 DO I ilo <> IF NDIM 0 DO p{{ J I }} DUP F@ p{{ ilo I }} F@ F+ F2/ FDUP psum{ I } F! F! LOOP psum{ func y{ I } F! THEN LOOP NDIM nfunc + TO nfunc ; Public: \ Multidimensional minimization of the function func(x) where x[0..ndim-1] \ is a vector in ndim dimensions, by the downhill simplex method of \ Nelder and Mead. Set the value NDIM before calling this function. The matrix \ p[0..ndim][0..ndim-1] is input. Its ndim+1 rows are ndim-dimensional vectors \ which are the vertices of the starting simplex. Also input is the vector \ y[0..ndim], whose components must be pre-initialized to the values of func \ evaluated at the ndim+1 vertices (rows) of p; and ftol the fractional \ convergence tolerance to be achieved in the function value (*). \ On output, p and y will have been reset to ndim+1 new points all within ftol \ of a minimum function value, and nfunc gives the number of function \ evaluations taken (nfunc is made -1 to signal non-convergence). \ (*) NB fractional: When the minimum is 0.0 you have a problem! : amoeba ( 'p 'y xt -- steps ) ( F: tol -- ) FRAME| a | ( a=ftol ) defines func & y{ &! & p{{ &! & psum{ NDIM }malloc malloc-fail? & ptrY{ NDIM }malloc malloc-fail? OR ABORT" amoeba :: out of memory" 0 TO nfunc geomean BEGIN TEST-HOOK order'm fract-range a F< IF & ptrY{ }free & psum{ }free |FRAME nfunc EXIT THEN nfunc nfuncmax > IF & ptrY{ }free & psum{ }free |FRAME -1 EXIT THEN reflect FDUP y{ ilo } F@ F<= IF FDROP double FDROP ELSE ( ytry) y{ inhi } F@ F>= IF y{ ihi } F@ halve F<= IF store.p geomean THEN THEN THEN AGAIN ; Reset_Search_Order TEST-CODE? [IF] \ ------------------------------------------------------------ FLOAT DMATRIX p{{ \ ndim+1, ndim : ndim+1 n-dimensional vectors FLOAT DMATRIX y{ \ the ndim+1 results when func evaluates the p vectors \ The function we'll seek the minimum of. : PARABOLOID ( 'mat -- ) ( F: -- r ) LOCALS| m{ | 1e NDIM 0 DO m{ I } F@ FSQR F+ LOOP ; \ Read ahead in a text file. This doesn't work with a terminal. \ A nice feature: the read text is interpreted, so { 1 2 + } works! \ Data starts on the _next_ line. : READ-INFILE REFILL 0= ABORT" REFILL :: not possible" SOURCE EVALUATE ; \ Read a matrix (won't work from the terminal). : }}FREAD ( 'head rows cols -- ) LOCALS| cols rows m{{ | m{{ rows cols }}malloc malloc-fail? ABORT" }}FREAD failed" m{{ EXECUTE TO m{{ rows 0 ?DO READ-INFILE ( coefficients) 0 cols 1- DO m{{ J I }} F! -1 +LOOP LOOP REFILL 0= ABORT" REFILL :: not possible" ; 3 TO NDIM FLOAT DMATRIX pp{{ & pp{{ NDIM 1+ NDIM }}FREAD \ array 0..ndim, 0..ndim-1 10e 0e 0e 0e 10e 0e 0e 0e 10e 10e 10e 10e FLOAT DARRAY m'{ : INIT-MATRICES & p{{ NDIM 1+ NDIM }}malloc pp{{ p{{ NDIM 1+ NDIM }}fcopy & m'{ NDIM }malloc & y{ NDIM 1+ }malloc NDIM 1+ 0 DO NDIM 0 ?DO p{{ J I }} F@ m'{ I } F! LOOP m'{ PARABOLOID y{ I } F! LOOP & m'{ }free ; : CLEANUP ( -- ) & y{ }free & p{{ }}free ; : TEST-AMOEBA ( F: tol -- ) INIT-MATRICES p{{ y{ use( PARABOLOID ( tol) FDUP amoeba CR . ." steps taken by AMOEBA when tol is " F. CR NDIM DUP 1+ . ." points in " 0 .R ." -dim space:" NDIM 1+ NDIM p{{ CR }}fprint CR ." Function values at these points :" NDIM 1+ y{ CR }fprint CLEANUP ; : SHOW-PROBLEM INIT-MATRICES CR NDIM DUP 1+ . ." points in " 0 .R ." -dim space:" NDIM 1+ NDIM p{{ CR }}fprint CR ." Function values at these points :" NDIM 1+ y{ CR }fprint CLEANUP ; : .ABOUT CR ." Type SHOW-PROBLEM for the initial problem." CR ." Try TEST-AMOEBA " CR CR ." is a fractional tolerance, take care when the minimum=0" ; [ELSE] : .ABOUT CR ." F: AMOEBA <#iters> " CR CR ." is a matrix of ndim+1 vectors with ndim components" CR ." is a matrix with the ndim+1 values of FUNC at these vectors" CR ." is the address of FUNC -- <> F: <> -- " CR ." where is an array of ndim variables" CR ." is a fractional tolerance, take care when the minimum=0" CR ." <#iters> is -1 when no convergence could be reached." CR CR ." Be sure to set NDIM before calling AMOEBA." ; [THEN] ( * End of Source * )