Arjen Markus (3 may 2005) I came across the idea of automatic differentiation in a Fortran 90 program by Simon Geard that he was writing based on a Python example for the Language Shootout.
Larry Smith This seems like an excellent application for 'C [https://www.researchgate.net/publication/2474569_C_and_tcc_A_Language_and_Compiler_for_Dynamic_Code_Generation] . This might make an excellent addition to Tcl. [L1 ] [L2 ] [L3 ]
While in Fortran 90 it is possible to define operators for so-called derived types, we can not do this in Tcl, at least not directly. The technique however is quite accessible and this page shows one way of doing it.
First let me explain what the idea is:
Symbolic manipulation is a kind of computer algebra that allows the computer to manipulate mathematical expressions, like rearranging the terms in the expression to make it more concise and easier to understand:
a**2 + 2*a*b + b**2 ==> (a+b)**2
One kind of manipulation is to find the derivative of some function:
f(x) = x**2 + 2*x df --(x) = 2*x + 2 dx
The rules for that kind of manipulation are simple and can applied mechanically (in contrast: finding the primitive is a challenge that quite often leads to completely new classes of functions!)
A very nice technique is to let the expression do the hard work: this is called automatic differentiation and here is one way to do it in Tcl:
proc f {x} {
return [+ [* $x $x] [* 2 $x]]
}set x [list 2.0 1.0]
So, the rules:
f(x) = g(x) + h(x) ==> f'(x) = g'(x) + h'(x) f(x) = g(x) * h(x) ==> f'(x) = g'(x)*h(x) + g(x)*h'(x)
etc. can be computed without us having to do more than specifying the way operators like + and * behave under differentiation.
namespace eval ::derivs {
proc + {x y} {
if { [llength $x] == 1 } {
set x [list $x 0.0]
}
if { [llength $y] == 1 } {
set y [list $y 0.0]
}
set v [expr {[lindex $x 0]+[lindex $y 0]}]
set dv [expr {[lindex $x 1]+[lindex $y 1]}]
return [list $v $dv]
}
proc * {x y} {
if { [llength $x] == 1 } {
set x [list $x 0.0]
}
if { [llength $y] == 1 } {
set y [list $y 0.0]
}
set v [expr {[lindex $x 0]*[lindex $y 0]}]
set dv [expr {[lindex $x 1]*[lindex $y 0] +
[lindex $x 0]*[lindex $y 1]}]
return [list $v $dv]
}
# - and / are left as an exercise ...
#
# Also define a few functions to make it more interesting
#
foreach {f df} {sin cos exp exp cos -sin log 1.0/ sqrt 0.5/sqrt} {
proc $f {x} [string map [list F $f DF $df] {
if { [llength $x] == 1 } {
set x [list $x 0.0]
}
set v [expr {F([lindex $x 0])}]
set dv [expr {[lindex $x 1]*DF([lindex $x 0])}]
return [list $v $dv]
}]
}
}
#
# A few tests:
# Define two (global) functions
#
proc f {x} {
return [* $x [+ 2 $x]]
}
proc g {x} {
return [* 2 [* [sin $x] [cos $x]]]
}
namespace eval ::derivs {
#
# f'(x) = (2+x)+x = 2+2x
# g'(x) = 2*cos(x)**2-2*sin(x)**2 = 2*cos(2x)
#
# "Import" the two functions - well, this is one way to do it :)
#
foreach p {f g} {
proc $p [info args ::$p] [info body ::$p]
}
#
# Test the derivatives
#
puts "Parabola:"
foreach x {0.0 0.5 1.0 1.5 2.0} {
set xl [list $x 1.0] ;# Very important!!
puts "[f $xl] -- [expr {2.0+2.0*$x}]"
}
puts "sin(2x):"
foreach x {0.0 0.5 1.0 1.5 2.0} {
set xl [list $x 1.0]
puts "[g $xl] -- [expr {2.0*cos(2.0*$x)}]"
}
}This can be extended to the evaluation of second derivatives if you like, or you can construct the expression that evaluates the derivative.
Another, somewhat simpler, variation is to determine the size of the expression tree - something that is important with genetic programming techniques. Merely evaluate the above functions in a different namespace ... a neat trick Donal Fellows mentioned on comp.lang.tcl:
namespace eval ::treesize {
foreach op {+ * - /} {
#
# Only the "arity" counts ...
#
proc $op {x y} {
set count 1
if { [llength $x] == 1 } {
incr count
} else {
incr count [lindex $x 1]
}
if { [llength $y] == 1 } {
incr count
} else {
incr count [lindex $y 1]
}
return [list COUNT $count]
}
}
#
# Also define a few functions to make it more interesting
#
foreach f {sin exp cos} {
proc $f {x} {
set count 1
if { [llength $x] == 1 } {
incr count
} else {
incr count [lindex $x 1]
}
return [list COUNT $count]
}
}
}
#
# Reuse the functions above ...
#
namespace eval ::treesize {
#
# "Import" the two functions
#
foreach p {f g} {
proc $p [info args ::$p] [info body ::$p]
}
puts "Size of 'f': [f 1]"
puts "Size of 'g': [g 1]"
}So, with the proper operators and mathematical functions we need only to define the actual expression once to automatically derive all kinds of mathematical and symbolic properties!
Suppose you have an expression like:
f(x) = ln( x + sqrt(1+x**2) )
and you need to compute the derivative, that is trivial with the code above. Manually verifying that its first derivative is:
1
f'(x) = ------------
sqrt(1+x**2)and then computing the value would take some time :)
TV From an earlier look at this page I recall the algorithm has limitation with respect to the types of forms of equations? If it is a general derivative routine, it would be good, regardless of the formula format.
Don't forget Maxima can also do derivations symbolically, quite powerfully, and even the reverse, which is a lot harder still, and it can be driven by tcl.
AM The method presented above is quite general, the implementation may be a bit limited due to my laziness :). The very attractive thing about it is that you do not need to resort to numerical differentiation which suffers from accuracy problems, not to mention problems with singularities ... This is an exact method with no caveats beyond the mathematical - I mean: the derivative has to exist.
SYMDIFF is a tool for symbolic differentiation. Features: