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

Mode: | factor |

Date: | Sun, 8 Feb 2009 02:56:32 |

USING: kernel math math.functions math.derivatives accessors words generalizations sequences generic.parser fry locals compiler.units continuations quotations combinators macros make ; IN: math.dual TUPLE: dual ordinary-part epsilon-part ; C: <dual> dual ! Ordinary numbers implement the dual protocol by returning themselves as the ! ordinary part, and 0 as the epsilon part. M: number ordinary-part>> ; M: number epsilon-part>> drop 0 ; : unpack-dual ( dual -- ordinary-part epsilon-part ) [ ordinary-part>> ] [ epsilon-part>> ] bi ; <PRIVATE ! :: chain-rule ( derivative-list n -- x ) ! { [ [ epsilon-part>> ] n napply ] [ [ ordinary-part>> ] n napply ] } ! n ncleave ! derivative-list [ n ncurry ] n nwith map ! spread n narray sum ; inline : my-join ( seq glue -- seq ) [ '[ [ _ % ] [ , ] interleave ] ] keep make ; : [derivative-at-point] ( derivative-list -- quot ) [ <reversed> ] [ length '[ _ nkeep ] ] bi my-join [ call ] compose ; ! MACRO: derivative-at-point ( xincs vals derivative-list -- yincs ) ! [derivative-at-point] ; :: chain-rule ( derivative-list n -- x ) [ [ epsilon-part>> ] n napply ] n nkeep [ ordinary-part>> ] n napply derivative-list [derivative-at-point] call n narray sum ; inline PRIVATE> ! dual-op is similar to execute, but promotes the word to operate on duals. ! Uses the "derivative" word-prop, which holds a list of quotations giving the ! partial derivative of the word with respect to each of its argumets. ! :: [dual-op] ( word -- quot ) ! word "derivative" word-prop :> derivative-list ! derivative-list length :> n ! [ [ [ ordinary-part>> ] n napply word execute ] n nkeep ! derivative-list n chain-rule ! <dual> ] ; :: [dual-op] ( word -- quot ) word "derivative" word-prop :> derivative-list derivative-list length :> n n word '[ [ ordinary-part>> ] _ napply _ execute ] n derivative-list n '[ _ _ nkeep _ _ chain-rule <dual> ] ; MACRO: dual-op ( word -- ) [dual-op] ; : define-dual-method ( word -- ) [ \ dual swap create-method-in ] keep [dual-op] define ; ! Specialize math functions to operate on dual numbers. << { sqrt exp log sin cos tan sinh cosh tanh atan } [ define-dual-method ] each >> ! Inverse methods { asin, acos, asinh, acosh, atanh } are not generic, ! so there is no way to specialize them for dual numbers. ! Arithmetic methods are not generic (yet?), so we have to define special ! versions of them to operate on dual numbers. : d+ ( x y -- x+y ) \ + dual-op ; : d- ( x y -- x+y ) \ - dual-op ; : d* ( x y -- x*y ) \ * dual-op ; : d/ ( x y -- x/y ) \ / dual-op ; : d^ ( x y -- x^y ) \ ^ dual-op ;

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

Mode: | factor |

Date: | Sun, 8 Feb 2009 03:07:24 |

! Copyright (C) 2009 Jason Merrill. ! See http://factorcode.org/license.txt for BSD license. USING: kernel math math.functions math.derivatives.parser ; IN: math.derivatives DERIVATIVE: + [ 2drop ] [ 2drop ] DERIVATIVE: - [ 2drop ] [ 2drop neg ] DERIVATIVE: * [ nip * ] [ drop * ] DERIVATIVE: / [ nip / ] [ sq / neg * ] ! Conditional checks if the epsilon-part of the exponent is 0 to avoid getting ! float answers for integer powers. DERIVATIVE: ^ [ [ 1 - ^ ] keep * * ] [ [ dup zero? ] 2dip [ 3drop 0 ] [ [ ^ ] keep log * * ] if ] DERIVATIVE: sqrt [ 2 * / ] DERIVATIVE: exp [ exp * ] DERIVATIVE: log [ / ] DERIVATIVE: sin [ cos * ] DERIVATIVE: cos [ sin neg * ] DERIVATIVE: tan [ sec sq * ] DERIVATIVE: sinh [ cosh * ] DERIVATIVE: cosh [ sinh * ] DERIVATIVE: tanh [ sech sq * ] DERIVATIVE: asin [ sq neg 1 + sqrt / ] DERIVATIVE: acos [ sq neg 1 + sqrt neg / ] DERIVATIVE: atan [ sq 1 + / ] DERIVATIVE: asinh [ sq 1 + sqrt / ] DERIVATIVE: acosh [ [ 1 + sqrt ] [ 1 - sqrt ] bi * / ] DERIVATIVE: atanh [ sq neg 1 + / ]

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

Mode: | factor |

Date: | Sun, 8 Feb 2009 03:08:06 |

! Copyright (C) 2009 Jason Merrill. ! See http://factorcode.org/license.txt for BSD license. USING: kernel parser words effects accessors sequences math.ranges ; IN: math.derivatives.parser : DERIVATIVE: scan-object dup stack-effect in>> length [1,b] [ drop scan-object ] map "derivative" set-word-prop ; parsing

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

Mode: | factor |

Date: | Sun, 8 Feb 2009 03:08:51 |

! Copyright (C) 2009 Jason Merrill. ! See http://factorcode.org/license.txt for BSD license. USING: kernel math math.functions math.dual assocs fry accessors ; IN: automatic-differentiation ! 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>> ] ; : 1diff ( x quot -- y ) [1diff] call ; inline

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

Mode: | factor |

Date: | Sun, 8 Feb 2009 03:09:38 |

! Copyright (C) 2009 Your name. ! See http://factorcode.org/license.txt for BSD license. USING: kernel tools.test automatic-differentiation math math.functions math.constants math.dual ; IN: automatic-differentiation.tests [ t ] [ 0 [ sin ] 1diff 1.0 epsilon 10 * ~ ] unit-test [ t ] [ 0 [ cos ] 1diff 0 epsilon 10 * ~ ] unit-test [ 1 ] [ 8 [ 3 + ] 1diff ] unit-test [ 3 5 ] [ 1 5 <dual> 2 d+ unpack-dual ] unit-test [ 0 -1 ] [ 1 5 <dual> 1 6 <dual> d- unpack-dual ] unit-test [ 2 1 ] [ 2 3 <dual> 1 -1 <dual> d* unpack-dual ] unit-test [ -1/4 ] [ 2 [ 1 swap / ] 1diff ] unit-test [ 2 ] [ 1 [ 2 ^ ] 1diff ] unit-test [ 1/8 ] [ 4 [ sqrt ] 1diff ] unit-test

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

Mode: | factor |

Date: | Sun, 8 Feb 2009 04:21:06 |

! Copyright (C) 2009 Jason Merrill. ! See http://factorcode.org/license.txt for BSD license. USING: kernel math math.functions math.derivatives accessors words generalizations sequences generic.parser fry locals compiler.units continuations quotations combinators macros make ; IN: math.dual TUPLE: dual ordinary-part epsilon-part ; C: <dual> dual ! Ordinary numbers implement the dual protocol by returning themselves as the ! ordinary part, and 0 as the epsilon part. M: number ordinary-part>> ; M: number epsilon-part>> drop 0 ; : unpack-dual ( dual -- ordinary-part epsilon-part ) [ ordinary-part>> ] [ epsilon-part>> ] bi ; <PRIVATE ! :: chain-rule ( derivative-list n -- x ) ! { [ [ epsilon-part>> ] n napply ] [ [ ordinary-part>> ] n napply ] } ! n ncleave ! derivative-list [ n ncurry ] n nwith map ! spread n narray sum ; inline MACRO: weave ( n -- quot ) [ '[ [ ] _ ncurry ] ] [ <reversed> [ '[ _ ndip ] ] map ] bi '[ @ _ cleave ] ; MACRO: nspread ( quots n -- quot ) over empty? [ 2drop [ ] ] [ [ [ but-last ] dip ] [ [ peek ] dip ] 2bi swap '[ [ _ _ nspread ] _ ndip @ ] ] if ; :: chain-rule ( derivative-list n -- x ) [ [ epsilon-part>> ] n napply ] n nkeep [ ordinary-part>> ] n napply n weave derivative-list n 1+ nspread n narray sum ; inline PRIVATE> ! dual-op is similar to execute, but promotes the word to operate on duals. ! Uses the "derivative" word-prop, which holds a list of quotations giving the ! partial derivative of the word with respect to each of its argumets. ! :: [dual-op] ( word -- quot ) ! word "derivative" word-prop :> derivative-list ! derivative-list length :> n ! [ [ [ ordinary-part>> ] n napply word execute ] n nkeep ! derivative-list n chain-rule ! <dual> ] ; :: [dual-op] ( word -- quot ) word "derivative" word-prop :> derivative-list derivative-list length :> n n word '[ [ ordinary-part>> ] _ napply _ execute ] n derivative-list n '[ _ _ nkeep _ _ chain-rule <dual> ] ; MACRO: dual-op ( word -- ) [dual-op] ; : define-dual-method ( word -- ) [ \ dual swap create-method-in ] keep [dual-op] define ; ! Specialize math functions to operate on dual numbers. << { sqrt exp log sin cos tan sinh cosh tanh atan } [ define-dual-method ] each >> ! Inverse methods { asin, acos, asinh, acosh, atanh } are not generic, ! so there is no way to specialize them for dual numbers. ! Arithmetic methods are not generic (yet?), so we have to define special ! versions of them to operate on dual numbers. : d+ ( x y -- x+y ) \ + dual-op ; : d- ( x y -- x+y ) \ - dual-op ; : d* ( x y -- x*y ) \ * dual-op ; : d/ ( x y -- x/y ) \ / dual-op ; : d^ ( x y -- x^y ) \ ^ dual-op ;

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

Author: | |

Mode: | |

Body: | |