Re: Differentiable Forth

I organized a bit the work.
here are 3 programs:
    - dual_numbers.fs, it containes operations and elementary functions
of dual_numbers
    - autodiff.fs, it includes dual_numbers.fs, it defines der,
gradient, jacobian
and - ad_tests.fs, it includes autodiff.fs, it containes some examples:
                   1 func of 1 var ---> der
                   1 func of 3 var ---> gradient
                   2 func of 3 vars ---> jacobian
                   3 func of 3 vars ---> jacobian
                   4 func of 2 vars ---> jacobian
                   4 func of 3 vars ---> jacobian
Here, begins the code for dial_numbers.fs
: -frot frot frot ;
: f>dl   ( f: a -- a 0) 0e ; \ real to dual value
: f>dl_d ( f: a -- a 1) 1e ; \ real to dual value with respect to
differentiate
: dl>rl  ( f: a b -- a) fdrop ;
: dl>eps ( f: a b -- b) fnip ;
: dl. ( f: a b -- ) fswap f. ." + eps " f. ;
: dl(,). ( f: a b -- ) ." (" fswap f. ." , " f. ." )" ;
: dl+ ( f: a b c d -- a+c b+d)
    frot ( f: a c d b)
    f+   ( f: a c b+d)
    -frot ( f: b+d a c)
    f+   ( f: b+d a+c)
    fswap ( f: a+c b+d) ;
: dl- ( f: a b c d -- a-c b-d)
    frot ( f: a c d b)
    fswap f-   ( f: a c b-d)
    -frot ( f: b-d a c)
    f-   ( f: b-d a-c)
    fswap ( f: a-c b-d) ;
fvariable b*c
: dl* ( f: a b c d -- a*c  a*d+b*c)
    -frot ( f: a d b c)
    ftuck f* ( f: a d c b*c)
    b*c f! ( f: a d c)
    frot ( f: d c a)
    ftuck ( f: d a c a)
    f*    ( f: d a a*c)
    -frot  ( f: a*c d a)
    f* b*c f@ f+ ( f: a*c a*d+b*c)
;
: 1/dl ( f: a b -- 1/a -b*1/a^2)
    fswap 1/f ( f: b 1/a)
    ftuck     ( f: 1/a b 1/a)
    fdup f*   ( f: 1/a b 1/a^2)
    f* fnegate ( f: 1/a -b/a^2)
;
: dl/ ( f: a b c d -- a/c rest)
    1/dl dl*
;
: dl^f ( f: a b c -- a^c b*c)
    ftuck ( f: a c b c)
    f* -frot ( f: b*c a c)
    f** fswap ( f: a^c b*c)
;
: dlnegate ( f: a b -- c d)
     fnegate fswap fnegate fswap
;
Here, dual_numbers.fs finishes
--------------------------------
Here, autodiff.fs begins
include dual_numbers.fs
\ -----------
\ vector valued functions of several variables
: funcarray ( n -- )
    create dup , cells allot
    does> dup @
;
: funcarray! ( func address size i --)
     nip cells + cell+ !
;
: funcarray@ ( address size i -- func)
     nip cells + cell+ @
;
\ derivatives
\ derivative 1 func of 1 var
variable func
: der 1e ' func ! func @ execute dl>eps ;
\ gradient and jacobian
: dlarray ( n --) \ does part ( -- address size)
    create dup , 2* floats allot does> dup @  ;
: dlarray! ( address size i --) ( dl: dl --)
    nip 2* floats + cell+ dl! ;
: dlarray@ ( address size i --) ( dl: -- dl)
    nip 2* floats + cell+ dl@ ;
variable func
: (gradient) ( address size --)
    ' func !
    dup 0 do
      dup 0 do
        2dup i dlarray@
        i j = if
          0e 1e dl+
        then
      loop
      func @ execute dl>eps
    loop
    2drop
;
: ([gradient]) ( address size --)
    dup 0 do
      dup 0 do
        2dup i dlarray@
        i j = if
          0e 1e dl+
        then
      loop
      func @ execute dl>eps
    loop
    2drop
;
: (.gradient)
    ' func !
    dup 0 do
      dup 0 do
        2dup i dlarray@
        i j = if
          0e 1e dl+
        then
      loop
      func @ execute dl>eps
      cr f.
    loop
    2drop
 ;
: ([.gradient])
    dup 0 do
      dup 0 do
        2dup i dlarray@
        i j = if
          0e 1e dl+
        then
      loop
      func @ execute dl>eps
      f.
    loop
    2drop
 ;
\ gradient
: .gradient 0 funcarray@ func ! ([.gradient]) ;
: gradient 0 funcarray@ func ! ([gradient]) ;
\ jacobian
2variable (funcarray)
2variable (dlarray)
: .jacobian
   (funcarray) 2!
   (dlarray) 2!
   (funcarray) 2@ nip 0 do
      (funcarray) 2@ i funcarray@ func !
      (dlarray) 2@ ([.gradient])
      cr
   loop
;
: jacobian
   (funcarray) 2!
   (dlarray) 2!
   (funcarray) 2@ nip 0 do
      (funcarray) 2@ i funcarray@ func !
      (dlarray) 2@ ([gradient])
   loop
;
Here, autodiff.fs finishes
-----------------------------
Here, the ad_tests.fs begins
include autodiff.fs
\ 1 func of 1 var
: dl_f() ( dl: x -- y) dldup dlsin dl* ; \ f(x) = x*sin(x)
\ for x = 1, df/fx = ?
cr 1e der dl_f() f.
: df/dx ( f: x -- y) fdup fdup fcos f* fswap fsin f+ ;
cr 1e df/dx f.
cr cr
\ 3 variables x, y and z in the array xyz
3 dlarray xyz
\ x = 1, y = 2 and z = 3
1e 0e xyz 0 dlarray!
2e 0e xyz 1 dlarray!
3e 0e xyz 2 dlarray!
\ 1 func of 3 vars
: dl_f() ( dl: x y z -- x + y * z) dl* dl+ ;
cr xyz (gradient) dl_f() frot f. fswap f. f.
cr xyz (.gradient) dl_f()
cr
\ -----------
\ vector valued functions
\ 1 func of 3 vars
1 funcarray f()
' dl_f() f() 0 funcarray!
cr xyz f() .gradient
\ cr xyz f() gradient \ the gradient values are in fpstack ( f: df/dx
df/dy df/dz)
\ 2 func of 3 var
: dl_f1() ( dl: x y z -- x+y*z) dl* dl+ ;
: dl_f2() ( dl: x y z -- x*y+z) dlrot dlrot dl* dl+ ;
2 funcarray f2()
' dl_f1() f2() 0 funcarray!
' dl_f2() f2() 1 funcarray!
cr xyz f2() .jacobian
\ cr xyz f2() jacobian \ the jacobian values are in fpstack ( f: df1/dx
df1/dy df1/dz df2/dx df2/dy df3/dz )
\ -----------
\ 3 func of 3 var
: dl_f1() ( dl: x y z -- x+y*z) dl* dl+ ;
: dl_f2() ( dl: x y z -- x*y+z) dlrot dlrot dl* dl+ ;
: dl_f3() ( dl: x y z -- y+z*x) dlrot dl* dl+ ;
3 funcarray f3()
' dl_f1() f3() 0 funcarray!
' dl_f2() f3() 1 funcarray!
' dl_f3() f3() 2 funcarray!
cr xyz f3() .jacobian
\ cr xyz f3() jacobian \ the jacobian values are in fpstack ( f: df1/dx
df1/dy df1/dz ... df3/dx df3/dy df3/dz )
\ -----------
\ 4 func of 2 var
2 dlarray xy
5e f>dl xy 0 dlarray!
6e f>dl xy 1 dlarray!
: dl_f1() ( dl: x y -- x+2*y) 2e f>dl dl* dl+ ;
: dl_f2() ( dl: x y -- x*y) dl* ;
: dl_f3() ( dl: x y -- x^y) dl^ ;
: dl_f4() ( dl: x y -- y^x) dlswap dl^ ;
4 funcarray f4()
' dl_f1() f4() 0 funcarray!
' dl_f2() f4() 1 funcarray!
' dl_f3() f4() 2 funcarray!
' dl_f4() f4() 3 funcarray!
cr xy f4() .jacobian
\ cr xy f4() jacobian \ the jacobian values are in fpstack ( f: df1/dx
df1/dy df1/dz ... df4/dx df4/dy df4/dz )
\ 4 func of 3 var
: dl_f1() ( dl: x y z -- x+y*z) dl* dl+ ;
: dl_f2() ( dl: x y z -- x*y+z) dlrot dlrot dl* dl+ ;
: dl_f3() ( dl: x y z -- y+z*x) dlrot dl* dl+ ;
: dl_f4() ( dl: x y z -- x*y*z) dl* dl* ;
4 funcarray f4()
' dl_f1() f4() 0 funcarray!
' dl_f2() f4() 1 funcarray!
' dl_f3() f4() 2 funcarray!
' dl_f4() f4() 3 funcarray!
cr xyz f4() .jacobian
\ cr xyz f4() jacobian \ the jacobian values are in fpstack ( f: df1/dx
df1/dy df1/dz ... df4/dx df4/dy df4/dz )
cr cr .( done) cr cr
Here, ad_tests.fs finishes
Ahmed