Author: | Jason Merrill |
---|---|

Mode: | factor |

Date: | Thu, 22 Jan 2009 01:09:28 |

! Copyright (C) 2009 Jason Merrill. ! See http://factorcode.org/license.txt for BSD license. USING: kernel fry accessors math math.functions locals assocs ; IN: automatic-differentiation TUPLE: dual ordinary-part epsilon-part ; C: <dual> dual M: number ordinary-part>> ; M: number epsilon-part>> drop 0 ; : unpack-dual ( dual -- ordinary-part epsilon-part ) [ ordinary-part>> ] [ epsilon-part>> ] bi ; : d+ ( x y -- x+y ) [ [ ordinary-part>> ] bi@ + ] [ [ epsilon-part>> ] bi@ + ] 2bi <dual> ; : dneg ( x -- -x ) unpack-dual [ neg ] bi@ <dual> ; : d- ( x y -- x-y ) dneg d+ ; : d* ( x y -- x*y ) [ [ ordinary-part>> ] bi@ * ] [ [ [ epsilon-part>> ] [ ordinary-part>> ] bi* * ] [ [ ordinary-part>> ] [ epsilon-part>> ] bi* * ] 2bi + ] 2bi <dual> ; ! Not written in terms of recip in order to avoid loss of precision. : d/ ( x y -- x/y ) [ [ ordinary-part>> ] bi@ / ] [ [ [ epsilon-part>> ] [ ordinary-part>> ] bi* / ] [ [ ordinary-part>> ] [ unpack-dual ] bi* swap sq / * ] 2bi - ] 2bi <dual> ; :: (d^) ( xo xe yo ye -- x^y ) xo yo ^ xo yo 1 - ^ yo * xe * ye 0 = [ ] [ xo yo ^ yo log * ye * + ] if <dual> ; : d^ ( x y -- x^y ) [ unpack-dual ] bi@ (d^) ; M: dual sin [ ordinary-part>> sin ] [ unpack-dual swap cos * ] bi <dual> ; M: dual cos [ ordinary-part>> cos ] [ unpack-dual swap sin neg * ] bi <dual> ; M: dual exp [ ordinary-part>> exp ] [ unpack-dual swap exp * ] bi <dual> ; M: dual log [ ordinary-part>> log ] [ unpack-dual swap / ] bi <dual> ; ! 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 <dual> @ 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.

Author: | jwmerrill |
---|---|

Mode: | factor |

Date: | Thu, 22 Jan 2009 06:38:30 |

! 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> 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 <dual> ; inline : 2dual-op ( dual1 dual2 quot d1quot d2quot -- dual ) [ '[ [ ordinary-part>> ] bi@ @ ] ] tri@ [ [ epsilon-part>> ] bi@ ] 4 narray 2cleave 2array [ 2array ] dip v. <dual> ; 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 <dual> ; ! 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 <dual> @ 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.

Summary: | |
---|---|

Author: | |

Mode: | |

Body: | |