Some rationale behind this package presented after the source code.
# expr.tcl
# Symbolic manipulation of expressions.
set __expr_tests 0
if $__expr_tests {
lappend auto_path .
source testsupp.tcl
}
package require atypes 1.0 ;# See [Algebraic Types]
package provide expr 1.0
# Some set operations.
# I had no TclX under Mandrake Linux and it was easier to write
# than to search and then compile.
proc setunion {a b} {
array set r {}
foreach x $a {
set r($x) 1
}
foreach y $b {
set r($y) 1
}
return [array names r]
}
proc setintersect {a b} {
array set t {}
foreach x $a {
set t($x) 1
}
set r {}
foreach y $b {
if {[info exists t($y)]} {
lappend r $y
}
}
return $r
}
# returns all elements from $a that is not in $b.
proc setsub {a b} {
set r {}
array set t {}
foreach x $b {
set t($x) 1
}
foreach y $a {
if {![info exists t($y)]} {
lappend r $y
}
}
return $r
}
# --------------------------------------------------------------------------------------------
# AST definition.
atype BinOp BPlus BMinus BMul BDiv
atype Fun FSin FCos FExp FSqrt
atype Expr {EConst c} {EVar v} {ENeg expr} {EBin op a b} {EFun fun expr}
# --------------------------------------------------------------------------------------------
# Construction. Actually, just constructor abbreviation for convenience.
proc ezero {} { return [EConst 0.0] }
proc eone {} { return [EConst 1.0] }
proc evar {v} {
return [EVar $v]
}
proc econst {c} {
return [EConst $c]
}
proc _ebin {op a b} {
return [EBin $op $a $b]
}
# construct eplus, eminus, emul and ediv.
foreach op {plus minus mul div} {
proc e$op {a b} " _ebin B[string toupper $op 0 0] \$a \$b"
}
# construct function abbreviation.
foreach f {sin cos exp sqrt} {
proc e$f {a} "return \[EFun F[string toupper $f 0 0] \$a\]"
}
# And checking.
proc eisconst {e} {
match $e {
{EConst _} { return 1 }
_ { return 0}
}
}
proc eepsilon {} { return 0.00000001 }
proc eclose {a b} { return [expr abs($a - $b)<[eepsilon]] }
proc eiszero {e} {
match $e {
{EConst $x} { return [eclose $x 0] }
_ { return 0}
}
}
proc eisone {e} {
match $e {
{EConst $x} { return [eclose $x 1] }
_ { return 0}
}
}
# --------------------------------------------------------------------------------------------
# Simplification.
# This is not as short as Haskell version, but works as well.
proc esimp {e} {
match $e {
{EConst _} { return $e }
{EVar _} { return $e }
{EBin BPlus $a $b} {
set sa [esimp $a]
set sb [esimp $b]
if {[eiszero $sa]} {
return $sb
}
if {[eiszero $sb]} {
return $sa
}
if {[eisconst $sa] && [eisconst $sb]} {
return [EConst [expr [efromconst $sa]+[efromconst $sb]]]
}
return [EBin BPlus $sa $sb]
}
{EBin BMinus $a $b} {
set sa [esimp $a]
set sb [esimp $b]
if {[eiszero $sb]} {
return $sa
}
if {[eisconst $sa] && [eisconst $sb]} {
return [EConst [expr [efromconst $sa]-[efromconst $sb]]]
}
return [EBin BMinus $sa $sb]
}
{EBin BMul $a $b} {
set sa [esimp $a]
set sb [esimp $b]
if {[eiszero $sa]} {
return [ezero]
}
if {[eiszero $sb]} {
return [ezero]
}
if {[eisone $sa]} {
return $sb
}
if {[eisone $sb]} {
return $sa
}
if {[eisconst $sa] && [eisconst $sb]} {
return [EConst [expr [efromconst $sa]*[efromconst $sb]]]
}
return [EBin BMul $sa $sb]
}
{EBin BDiv $a $b} {
set sa [esimp $a]
set sb [esimp $b]
if {[eiszero $sa]} {
return [ezero]
}
if {[eisone $sb]} {
return $sa
}
if {[eisconst $sa] && [eisconst $sb]} {
return [EConst [expr [efromconst $sa]/[efromconst $sb]]]
}
return [EBin BMul $sa $sb]
}
{EFun $fun $e} {
set se [esimp $e]
if {[eisconst $se]} {
match $fun {
FSin {set f sin}
FCos {set f cos}
FExp {set f exp}
FSqrt { set f sqrt}
}
return [EConst [expr $f([efromconst $se])]]
}
return [EFun $fun $se]
}
}
}
# --------------------------------------------------------------------------------------------
# Differentiation.
proc _diffbin {op a b v} {
set da [_ediff $a $v]
set db [_ediff $b $v]
match $op {
BMul {
return [eplus [emul $da $b] [emul $a $db]]
}
BDiv {
set num [eminus [emul $da $b] [emul $a $db]]
set den [emul $b $b]
return [ediv $num $den]
}
_ {
# Linear operations (+-).
_ebin $op [_ediff $a $v] [_ediff $b $v]
}
}
}
proc _difffun {fun e v} {
match $fun {
FSin { set d [ecos $e] }
FCos { set d [emul [econst -1] [esin $e]] }
FExp { set d [eexp $e] }
FSqrt { set d [ediv [const 0.5] [esqrt e]] }
}
return [emul [_ediff $e $v] $d]
}
proc _ediff {e v} {
match $e {
{EConst _} { ezero }
{EVar $x} {
if {[string equal $x $v]} {
eone
} else {
ezero
}
}
{EBin $op $a $b} {
_diffbin $op $a $b $v
}
{EFun $fun $e} {
_difffun $fun $e $v
}
}
}
proc ediff {e v} { return [esimp [_ediff $e $v]] }
# --------------------------------------------------------------------------------------------
# Different information about expression.
# Return list-set of variables from expression.
proc evars {e} {
match $e {
{EVar $v} { return [list $v] }
{EFun _ $a} { return [evars $a] }
{EBin _ $a $b} {
return [setunion [evars $a] [evars $b]]
}
_ { return {} }
}
}
# --------------------------------------------------------------------------------------------
# Evaluation (maybe, partial).
proc efromconst {e} {
match $e {
{EConst $x} { return $x }
_ { error "efromconst expects EConst!" }
}
}
proc _eeval {e context} {
match $e {
{EFun _ _} { error "No evaluation fo EFun!" }
{EBin $op $a $b} {
set ea [_eeval $a $context]
set eb [_eeval $b $context]
if {[eisconst $ea] && [eisconst $eb]} {
match $op {
BPlus { set op + }
BMinus { set op - }
BMul { set op * }
BDiv { set op / }
}
return [econst [expr [efromconst $ea] $op [efromconst $eb]]]
} else {
return [EBin $op $ea $eb]
}
}
{EVar $v} {
array set vals $context
if {[info exists vals($v)]} {
return [econst $vals($v)]
} else {
return $e
}
}
_ { return $e }
}
}
# (partially) evaluate expression given context as list of pairs ?variable value?..
proc eeval {e context} {
return [esimp [_eeval $e $context]]
}
# --------------------------------------------------------------------------------------------
# Tests
if $__expr_tests {
test {set a [eplus [evar "a"] [evar "b"]]}
test {set b [eexp [evar "a"]]}
test {set c [emul $a $b]}
test {evars $c}
test {set d [ediff $c a]}
test {set evres [eeval $a {a 1}]}
test {set evres [eeval $a {a 1 b 10}]}
}Why did I (SZ) need it?
I tried to solve some problem from Computer Vision and found that I need to find minimum of sum of quite complex squares. That could be done using Newton's method [L1 ] for multiple variables (using inverse Jacobian [L2 ] and all that). The extra problem was that formulation of those squares could change just when you turn around.
So I decided not to hardcode this quite complex sum of squares and its first derivative and Jacobian, but write this little package and differentiate on the fly.
We'll see how it turns out in practice.
BTW, I adapted this package from my set of expression operations written in Haskell.
Lars H: This appears to do symbolic differentiation. An alternative approach for your problem could be that of the Automatic differentiation page, where the function and its derivatives are all evaluated simultaneously, but no symbolic expression for the derivatives is made explicit.