! Copyright (C) 2009 Jason Merrill.
! See http://factorcode.org/license.txt for BSD license.
USING: kernel fry accessors math math.functions math.vectors
locals assocs arrays generalizations combinators ;
IN: automatic-differentiation
TUPLE: dual ordinary-part epsilon-part ;
C: dual
M: number ordinary-part>> ;
M: number epsilon-part>> drop 0 ;
: unpack-dual ( dual -- ordinary-part epsilon-part )
[ ordinary-part>> ] [ epsilon-part>> ] bi ;
! Takes a dual number, a quote that maps reals to reals, and another quote
! giving its derivative, and performs the analagous operation on the dual
! number.
: dual-op ( dual quot dquot -- dual' )
[ '[ ordinary-part>> @ ] ] [ '[ unpack-dual swap @ * ] ] bi*
bi ; inline
: 2dual-op ( dual1 dual2 quot d1quot d2quot -- dual )
[ '[ [ ordinary-part>> ] bi@ @ ] ] tri@
[ [ epsilon-part>> ] bi@ ]
4 narray 2cleave 2array [ 2array ] dip v. ; inline
! Is it worth generalizing the above to operations with arbitrary stack
! effect?
: d+ ( x y -- x+y ) [ + ] [ 2drop 1 ] [ 2drop 1 ] 2dual-op ;
: d- ( x y -- x-y ) [ - ] [ 2drop 1 ] [ 2drop -1 ] 2dual-op ;
: d* ( x y -- x*y ) [ * ] [ nip ] [ drop ] 2dual-op ;
! Careful about loss of precision with recip and *
: d/ ( x y -- x/y ) [ / ] [ nip recip ] [ sq / neg ] 2dual-op ;
! Gives floating point answers even for integer powers
: d^ ( x y -- x^y ) [ ^ ] [ [ 1 - ^ ] keep * ] [ [ ^ ] keep log * ] 2dual-op ;
! Given the above, is it worth allowing a custom way of combining values with
! increments? Mathematically speaking, it's always multiplication,
M: dual sin [ sin ] [ cos ] dual-op ;
M: dual cos [ cos ] [ sin neg ] dual-op ;
M: dual exp [ exp ] [ exp ] dual-op ;
! Beware loss of precision using recip followed by multiplication. Compare to
! explicit definition below.
M: dual log [ log ] [ recip ] dual-op ;
! M: dual log
! [ ordinary-part>> log ] [ unpack-dual swap / ] bi ;
! Replaces instances of math operators with their dual counterparts. This
! could be much smarter (work with nested quotations, inline word defs, etc.).
! If math operators could be overloaded like other generic words, there would
! be no need for this.
: lift ( quot -- quot' )
H{
{ + d+ }
{ - d- }
{ * d* }
{ / d/ }
{ ^ d^ }
}
substitute ;
: [1diff] ( quot -- quot' )
lift '[ 1 @ epsilon-part>> ] ; inline
: 1diff ( x quot -- y ) [1diff] call ; inline
! Plan for improvement:
! * Write a few unit tests.
! * Create a macro (?) define-derivative that takes a word and its derivative
! and specializes the word to operate on the dual class.
! * Similar to above for math-derivative.
! * Create a way to flatten definitions so that the math operators can be
! replaced with their dual counterparts.