Posted to tcl by auriocus at Thu May 15 19:19:06 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+$dy*$x}]
  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. puts "Argument: <$e>"
  56. lassign [dissectargs $e] x dx
  57. mkAD [expr {sin($x)}] [expr {cos($x)*$dx}]
  58. }
  59.  
  60. proc cos {e} {
  61. lassign [dissectargs $e] x dx
  62. mkAD [expr {cos($x)}] [expr {-sin($x)*$dx}]
  63. }
  64.  
  65. }
  66.  
  67. # now inject into VecTcl by renaming the numarray procedures
  68. namespace eval savedops {}
  69.  
  70. foreach op {+ - * / exp sin cos} {
  71. # this example throws them away....
  72. rename numarray::$op savedops::$op
  73. interp alias {} numarray::$op {} ADeval::$op
  74. }
  75.  
  76. # test
  77. # dependent variable x
  78. #
  79. set rx [expr {3.1415926535/6}] ;# 30°
  80. set x [mkAD $rx 1.0]
  81.  
  82. puts [vexpr {sin(x)}]
  83. puts "f=[expr {sin($rx)}], df/dx = [expr {cos($rx)}]"
  84.  
  85. # arjens example
  86. set a -0.5
  87. foreach rx {0 1 2 3 4 5} {
  88. set x [mkAD $rx 1.0]
  89.  
  90. puts [vexpr {exp(a*x)}]
  91. puts "f=[expr {exp($a*$rx)}] df=[expr {$a*exp($a*$rx)}]"
  92. }
  93.