Paste: dual numbers

Author:	Jason Merrill
Mode:	factor
Date:	Sun, 8 Feb 2009 02:56:32

Plain Text |

Delete Paste

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 ;

Annotation: math.derivatives

Author:	jwmerrill
Mode:	factor
Date:	Sun, 8 Feb 2009 03:07:24

Plain Text |

Delete Annotation

! 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 + / ]

Annotation: derivative-parser

Author:	jwmerrill
Mode:	factor
Date:	Sun, 8 Feb 2009 03:08:06

Plain Text |

Delete Annotation

! 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

Annotation: automatic-differentiation

Author:	jwmerrill
Mode:	factor
Date:	Sun, 8 Feb 2009 03:08:51

Plain Text |

Delete Annotation

! 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

Annotation: automatic-differentiation-tests

Author:	jwmerrill
Mode:	factor
Date:	Sun, 8 Feb 2009 03:09:38

Plain Text |

Delete Annotation

! 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

Annotation: dual numbers

Author:	jwmerrill
Mode:	factor
Date:	Sun, 8 Feb 2009 04:21:06

Plain Text |

Delete Annotation

! 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 ;

New Annotation

Summary:
Author:
Mode:	textPostgreSQLPowerCenter Parameter Fileactionscriptadaada95ans-forth94antantlrantlr4apacheconfapdlapplescriptaspaspect-jassembly-agcassembly-agsassembly-arm32assembly-m68kassembly-macro32assembly-mcs51assembly-parrotassembly-r2000assembly-x86avroawkbbatchbbjbcelbeanshellbibtexbinsource-agccc#c++cfscriptchillcilclipsclojurecmakecobolcoffeescriptcoldfusioncplex-lpcsscsvcvs-commitddartdjangodockerfiledotdoxygendrawj2ddsssleiffelembperlerlangfactorfhtmlforthfortranfortran90foxprofreemarkergcbasicgettextgherkingnuplotgogradlegroovygsphaskellhaxehexhl7v2hlslhtaccesshtmlhxmli4glicalendariconidlinforminiinno-setupinterlisiojamonjavajava module-infojavaccjavafxjavascriptjcljedit-actionsjflexjhtmljmkjsonjspkotlinlatexlexlilypondlispliterate-haskelllogslogtalklotosluam4macroschedulermailmakefilemaplemarkdownmavenmavscriptmetapostmlmodula3moinmqscmustache-templatemxmlmyghtymysqln3netrexxnqcnsis2objective-cobjectrexxoccamomnimarkopensipsoutlinepascalpatchperlphppikepl-sqlpl-sql9pl1plaintexplantumlpop11postscriptpovraypowerdynamopowershellprogressprologpropertiesprotobufpspptlpurepvwavepyrexpythonqdocrakefilercprdrebolredcoderelax-ng-compactrenderman-ribrestrfcrhtmlroffrpm-specrtfrubyrustrviews#s+sassbtscalaschemesdl/prsgmlshellscriptshtmlsipsippslateslaxsmalltalksmartysmi-mibsparqlsql-loadersqrsquidconfstatasvn-commitswiftswigtcltemplate-toolkittextexinfotldtransact-sqltsptwikityposcripturluscriptvalavbscriptvelocityverilogvhdlvisualbasicvrml2wellknowntextxmlxqxslyabyamlzpt
Body:

Done