Paste: automatic-differentiation

Author:	Jason Merrill
Mode:	factor
Date:	Thu, 22 Jan 2009 01:09:28

Plain Text |

Delete Paste

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

Annotation: factored out dual opreation definition

Author:	jwmerrill
Mode:	factor
Date:	Thu, 22 Jan 2009 06:38:30

Plain Text |

Delete Annotation

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

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