Posted to tcl by auriocus at Thu May 15 19:55:39 GMT 2014view raw

  1. # simple example for using vexpr to do automatic differentation (AD)
  2. # AD is a technique to compute an expression alongside with it's
  3. # analytic derivative using the chainrule
  4. package require vectcl
  5. namespace import vectcl::*
  6.  
  7. proc isAD {x} { expr {[lindex $x 0] eq "AD"} }
  8.  
  9. proc mkAD {x dx} { list AD $x $dx }
  10.  
  11. namespace eval ADeval {
  12. # the definition of an AD expression is
  13. # a list with AD, then x, then dx (only scalar for clarity)
  14. # also there is no forwarding to numarray
  15. proc dissectargs {args} {
  16. # result is pairwise, value and derivative
  17. set result {}
  18. foreach el $args {
  19. if {[isAD $el]} {
  20. lassign $el AD x dx
  21. lappend result $x $dx
  22. } else {
  23. lappend result $el 0
  24. }
  25. }
  26. return $result
  27. }
  28.  
  29. proc + {e1 e2} {
  30. lassign [dissectargs $e1 $e2] x dx y dy
  31. mkAD [expr {$x+$y}] [expr {$dx+$dy}]
  32. }
  33.  
  34. proc - {e1 e2} {
  35. lassign [dissectargs $e1 $e2] x dx y dy
  36. mkAD [expr {$x-$y}] [expr {$dx-$dy}]
  37. }
  38.  
  39. proc / {e1 e2} {
  40. lassign [dissectargs $e1 $e2] x dx y dy
  41. mkAD [expr {$x/$y}] [expr {($dx*$y - $dy*$x)/($y*$y)}]
  42. }
  43.  
  44. proc * {e1 e2} {
  45. lassign [dissectargs $e1 $e2] x dx y dy
  46. mkAD [expr {$x*$y}] [expr {$x*$dy+$dx*$y}]
  47. }
  48.  
  49. proc exp {e} {
  50. lassign [dissectargs $e] x dx
  51. mkAD [expr {exp($x)}] [expr {exp($x)*$dx}]
  52. }
  53.  
  54. proc sin {e} {
  55. lassign [dissectargs $e] x dx
  56. mkAD [expr {sin($x)}] [expr {cos($x)*$dx}]
  57. }
  58.  
  59. proc cos {e} {
  60. lassign [dissectargs $e] x dx
  61. mkAD [expr {cos($x)}] [expr {-sin($x)*$dx}]
  62. }
  63.  
  64. }
  65.  
  66. # now inject into VecTcl by renaming the numarray procedures
  67. namespace eval savedops {}
  68.  
  69. foreach op {+ - * / exp sin cos} {
  70. # this example throws away the numarray commands
  71. rename numarray::$op savedops::$op
  72. interp alias {} numarray::$op {} ADeval::$op
  73. }
  74.  
  75. # test
  76. # dependent variable x
  77. set rx [expr {3.1415926535/6}] ;# 30°
  78. set x [mkAD $rx 1.0]
  79.  
  80. puts [vexpr {sin(x)}]
  81. puts "f=[expr {sin($rx)}], df/dx = [expr {cos($rx)}]"
  82.  
  83. # Arjens example
  84. set a -0.5
  85. foreach rx {0 1 2 3 4 5} {
  86. set x [mkAD $rx 1.0]
  87. puts "x=$rx"
  88. puts "Automatic: [vexpr {exp(a*x)}]"
  89. puts "Manual: [mkAD [expr {exp($a*$rx)}] [expr {$a*exp($a*$rx)}]]"
  90.  
  91. puts "Automatic: [vexpr {2*sin(x)/cos(x)}]"
  92. puts "Manual: [mkAD [expr {2*sin($rx)/cos($rx)}] [expr {2.0/cos($rx)**2}]]"
  93. }