\ Simulated annealing minimizer
\ version of September 28th, 2000 - 20:23
\   2-dimensional case

FALSE [IF]
    Algorithm:
        0. Set the "temperature" T.
        1. Choose a point at random and compute f(x)
        2. If f < f_old then keep the new point;
        3. If f > f_old then compute P = exp[- (f - f_old)/T]
        4. Compute a random number -- if r < P then keep the new
           point else discard it.
        5. Reduce T gradually according to a reasonable schedule
           until the system settles down.

[THEN]

marker -anneal

: [undefined]   BL WORD  FIND  NIP  0=  ;

[undefined] f"     [IF]  include  ftran111.f  [THEN]
[undefined]  prng  [IF]  include  prng.f      [THEN]

v: fdummy

100 VALUE  Ntimes

: fvariables     0 DO  FVARIABLE  LOOP  ;

7  fvariables  T dT  oldF  xp yp  oldX  oldY

: transform    ( f: xi -- log[xi/[1-xi]] )
    1e0   FOVER  F-  F/  FLN  ;

: newXY ( f: -- x y)      \ x, y in (-inf to +inf)
    prng  transform  FDUP  xp F!
    prng  transform  FDUP  yp F!
;

0.1  seed 2!

: initialize   ( xt --)    ( f: T --)
    defines fdummy
    T  F!
    f" dT = T/10"
    newXY  fdummy  oldF F!
;

: keep   ( f: f[x,y] --)
    oldF F!   f" oldX = xp"   f" oldY = yp"
;

: select    ( f: f -- )
    FDUP                    ( f: f f)
    oldF F@  F-   FDUP        ( f: f df df)
    F0< IF    FDROP  keep   \ f < oldF so store f
              EXIT
    THEN

    T F@  F/  FNEGATE  FEXP ( f: f P)   \ f > oldF so compute P
    prng  F>  IF    keep                \ store f with prob P
              ELSE  FDROP  THEN         \ otherwise do nothing
;


: )anneal   ( xt --)    ( f: T --)
    initialize          ( --) ( f: --)
    BEGIN   T F@   F0>
    WHILE   Ntimes 0 DO
                newXY  fdummy  select
            LOOP
            f" T = T - dT"
    REPEAT              ( f: )
                        \ output
    f" oldF"   f" oldY"  f" oldX"
    CR ." (x,y) = " F. F.
    CR ." f_min = " F.
;

3 fvariables  r x y

: HW    ( f: x y -- f)
    y F!  x F!
    f" r = sqrt(x^2+y^2)"
    r F@  0.1e  F<
    IF  990e0  ELSE  f" 1/r^3 - 1/r + r + x"  THEN
;