unmaintained/math/derivatives/derivatives.factor

   1 ! Copyright (c) 2008 Reginald Keith Ford II, Eduardo Cavazos.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: kernel continuations combinators sequences math math.order math.ranges
   4     accessors float-arrays ;
   5 IN: math.derivatives
   6
   7 TUPLE: state x func h err i j errt fac hh ans a done ;
   8
   9 : largest-float ( -- x ) HEX: 7fefffffffffffff bits>double ; foldable
  10 : ntab ( -- val ) 8 ; inline
  11 : con ( -- val ) 1.6 ; inline
  12 : con2 ( -- val ) con con * ; inline
  13 : big ( -- val ) largest-float ; inline
  14 : safe ( -- val ) 2.0 ; inline
  15
  16 ! Yes, this was ported from C code.
  17 : a[i][i]     ( state -- elt ) [ i>>     ] [ i>>     ] [ a>> ] tri nth nth ;
  18 : a[j][i]     ( state -- elt ) [ i>>     ] [ j>>     ] [ a>> ] tri nth nth ;
  19 : a[j-1][i]   ( state -- elt ) [ i>>     ] [ j>> 1 - ] [ a>> ] tri nth nth ;
  20 : a[j-1][i-1] ( state -- elt ) [ i>> 1 - ] [ j>> 1 - ] [ a>> ] tri nth nth ;
  21 : a[i-1][i-1] ( state -- elt ) [ i>> 1 - ] [ i>> 1 - ] [ a>> ] tri nth nth ;
  22
  23 : check-h ( state -- state )
  24     dup h>> 0 = [ "h must be nonzero in dfridr" throw ] when ;
  25
  26 : init-a     ( state -- state ) ntab [ ntab <float-array> ] replicate >>a ;
  27 : init-hh    ( state -- state ) dup h>> >>hh ;
  28 : init-err   ( state -- state ) big >>err ;
  29 : update-hh  ( state -- state ) dup hh>> con / >>hh ;
  30 : reset-fac  ( state -- state ) con2 >>fac ;
  31 : update-fac ( state -- state ) dup fac>> con2 * >>fac ;
  32
  33 ! If error is decreased, save the improved answer
  34 : error-decreased? ( state -- state ? ) [ ] [ errt>> ] [ err>> ] tri <= ;
  35
  36 : save-improved-answer ( state -- state )
  37     dup err>>   >>errt
  38     dup a[j][i] >>ans ;
  39
  40 ! If higher order is worse by a significant factor SAFE, then quit early.
  41 : check-safe ( state -- state )
  42     dup [ [ a[i][i] ] [ a[i-1][i-1] ] bi - abs ]
  43     [ err>> safe * ] bi >= [ t >>done ] when ;
  44
  45 : x+hh ( state -- val ) [ x>> ] [ hh>> ] bi + ;
  46 : x-hh ( state -- val ) [ x>> ] [ hh>> ] bi - ;
  47
  48 : limit-approx ( state -- val )
  49     [
  50         [ [ x+hh ] [ func>> ] bi call ]
  51         [ [ x-hh ] [ func>> ] bi call ] bi -
  52     ] [ hh>> 2.0 * ] bi / ;
  53
  54 : a[0][0]! ( state -- state )
  55     { [ ] [ limit-approx ] [ drop 0 ] [ drop 0 ] [ a>> ] } cleave nth set-nth ;
  56
  57 : a[0][i]! ( state -- state )
  58     { [ ] [ limit-approx ] [ i>> ] [ drop 0 ] [ a>> ] } cleave nth set-nth ;
  59
  60 : a[j-1][i]*fac ( state -- val ) [ a[j-1][i] ] [ fac>> ] bi * ;
  61
  62 : new-a[j][i] ( state -- val )
  63     [ [ a[j-1][i]*fac ] [ a[j-1][i-1] ] bi - ]
  64     [ fac>> 1.0 - ] bi / ;
  65
  66 : a[j][i]! ( state -- state )
  67     { [ ] [ new-a[j][i] ] [ i>> ] [ j>> ] [ a>> ] } cleave nth set-nth ;
  68
  69 : update-errt ( state -- state )
  70     dup [ [ a[j][i] ] [ a[j-1][i] ] bi - abs ]
  71     [ [ a[j][i] ] [ a[j-1][i-1] ] bi - abs ] bi max >>errt ;
  72
  73 : not-done? ( state -- state ? ) dup done>> not ;
  74
  75 : derive ( state -- state )
  76     init-a
  77     check-h
  78     init-hh
  79     a[0][0]!
  80     init-err
  81     1 ntab [a,b) [
  82         >>i not-done? [
  83             update-hh
  84             a[0][i]!
  85             reset-fac
  86             1 over i>> [a,b] [
  87                 >>j
  88                 a[j][i]!
  89                 update-fac
  90                 update-errt
  91                 error-decreased? [ save-improved-answer ] when
  92             ] each check-safe
  93         ] when
  94    ] each ;
  95
  96 : derivative-state ( x func h err -- state )
  97     state new
  98     swap >>err
  99     swap >>h
 100     swap >>func
 101     swap >>x ;
 102
 103 ! For scientists:
 104 ! h should be .001 to .5 -- too small can cause bad convergence,
 105 ! h should be small enough to give the correct sgn(f'(x))
 106 ! err is the max tolerance of gain in error for a single iteration-
 107 : (derivative) ( x func h err -- ans error )
 108     derivative-state derive [ ans>> ] [ errt>> ] bi ;
 109
 110 : derivative ( x func -- m ) 0.01 2.0 (derivative) drop ;
 111 : derivative-func ( func -- der ) [ derivative ] curry ;