# An implementation of a a Bayesian belief network by way of
 # support-except and evidence-except algorithm by Jason Tang
 # ([email protected]).  The network is limited to just polytrees.
 
 # Populates the table of calculated beliefs using values given by the
 # Bayesian network.
 proc initBeliefs {} {
     array unset ::beliefs
     foreach n $::bayes(nodes) {
         if {$::bayes($n:p) == {}} {
             # head node
             set ::beliefs(+$n) $::bayes($n:w)
             set ::beliefs(~$n) [not $::bayes($n:w)]
         } else {
             # iterate through the parents and populate beliefs with
             # given evidence
             set weightNum 0
             foreach given [permutate $::bayes($n:p)] {
                 set given [join $given {}]
                 set ::beliefs(+$n|$given) [lindex $::bayes($n:w) $weightNum]
                 set ::beliefs(~$n|$given) [not $::beliefs(+$n|$given)]
                 incr weightNum
             }
         }
     }
 }
 
 # Given a child node and some condition by the parents, calculates the
 # probability of that child existing.  $a must be of form "+X" (to
 # indicate that state X should be true) or "~X" (to indicate
 # falsehood).  $given may be any complex truthhood (e.g., "+X~Y~Z") as
 # long as all of the states are the same distance away from $a.
 proc support {a given} {
     if {$given == {}} {
         # asked for a conditional variable, but a has no parents
         return {}
     }
     if {![info exists ::beliefs($a|$given)]} {
         set x [string range $a 1 end]
         set tmp $given
         while {[regexp {\A.([^+~]+)(.*)} $tmp foo g tmp]} {
             lappend gList $g
         }
         if [elemSubset $gList $::bayes($x:p)] {
             # condition x for just $given and no other variables
             # within x's parentage
             set prob 0
             foreach xgiven [permutate $::bayes($x:p)] {
                 set xgiven [join $xgiven {}]
                 if [truthSubset $given $xgiven] {
                     set p $::beliefs(+$x|$xgiven)
                     foreach r [truthRemainder $given $xgiven] {
                         set p [expr {double($p * [evidence $r])}]
                     }
                     set prob [expr {$prob + $p}]
                 }
             }
         } else {
             # search through parents until a path is found
             set foundPath 0
             foreach p $::bayes($x:p) {
                 set prob0 [support +$p $given]
                 if {$prob0 != {}} {
                     set prob1 [support ~$p $given]
                     set prob [expr {[support +$x +$p] * $prob0 +
                                     [support +$x ~$p] * $prob1}]
                     set foundPath 1
                     break
                 }
             }
             if !$foundPath {
                 error "no path found from $x to $given"
             }            
         }
         set ::beliefs(+$x|$given) $prob
         set ::beliefs(~$x|$given) [not $prob]
     }
 #    puts "support returns for $a|$given: $::beliefs($a|$given)"
     return $::beliefs($a|$given)
 }
 
 # Calculates the probability the truthhood of some state.  $a must be
 # of form "+X" (to indicate that state X should be true) or "~X" (to
 # indicate falsehood).
 proc evidence {a} {
     if {![info exists ::beliefs($a)]} {
         set x [string range $a 1 end]
         set prob 0.0
         foreach given [permutate $::bayes($x:p)] {
             set p [support +$x [join $given {}]]
             foreach g $given {
                 set p [expr {double($p * [evidence $g])}]
             }
             set prob [expr {$prob + $p}]
         }
         set ::beliefs(+$x) $prob
         set ::beliefs(~$x) [not $prob]
     }
 #    puts "evidence returns for $a: $::beliefs($a)"
     return $::beliefs($a)
 }
 
 # Uses Bayes' theorem to infer some parent condition ($a) given the
 # truthhood of a single child state ($given).  Both $a and $given must
 # be a single state of form "+X" (to indicate that state X should be
 # true) or "~X" (to indicate falsehood).
 proc infer {a given} {
     return [expr {[support $given $a] * [evidence $a] / [evidence $given]}]
 }
 
 ######################################################################
 # Bayesian utilities
 
 # Given a list of nodes, returns a permutation of them alternating
 # between true and false.  For example, given the list "a b c" this
 # function returns:
 #   +a+b+c ~a+b+c +a~b+c ~a~b+c +a+b~c ~a+b~c +a~b~c ~a~b~c
 # meaning {a, b, c}, {not a, b, c}, {a, not b, c}, etc.
 proc permutate {nodes} {
     if {$nodes == {}} {
         error "permutated across nothing!"
     }
     if {[llength $nodes] == 1} {
         return [list [list "+$nodes"] [list "~$nodes"]]
     } else {
         set head [lindex $nodes 0]
         foreach n [permutate [lrange $nodes 1 end]] {
             lappend retlist [concat "+$head" $n] [concat "~$head" $n]
         }
         return $retlist
     }
 }
 
 # returns non-zero if every element in needle is in the haystack list
 proc elemSubset {needle haystack} {
     foreach n $needle {
         if {[lsearch -exact $haystack $n] == -1} {
             return 0
         }
     }
     return 1
 }
 
 # Returns true if the truth $needle (e.g., "+a~b") is a subset of the
 # truth $haystack (e.g., "+a~x+y+b~z")
 proc truthSubset {needle haystack} {
     while {[regexp {\A(.[^+~]+)(.*)} $haystack foo h haystack]} {
         lappend stackList $h
     }
     while {[regexp {\A(.[^+~]+)(.*)} $needle foo n needle]} {
         if {[lsearch -exact $stackList $n] == -1} {
             return 0
         }
     }
     return 1
 }
 
 # Given a truth $needle (e.g., "+a~b") and a truth $haystack (e.g.,
 # "+a~x+y+b~z") return the truth without any from $needle in list form
 # (in this case, {~x +y ~z})
 proc truthRemainder {needle haystack} {
     while {[regexp {\A(.[^+~]+)(.*)} $needle foo n needle]} {
         lappend needleList $n
     }
     set remainder ""
     while {[regexp {\A(.[^+~]+)(.*)} $haystack foo h haystack]} {
         if {[lsearch -exact $needleList $h] == -1} {
             lappend remainder $h
         }
     }
     return $remainder
 }
 
 # Takes a probability and returns its inverse.
 proc not {p} {
     return [expr {1.0 - $p}]
 }
 
 # adds a node to the network.  be aware that weights must be phrased a
 # certain way
 proc addNode {shortName parents children weights {fullName ""}} {
     if {$fullName == ""} {
         set fullName $shortName
     }
     set ::bayes($shortName) $fullName
     set ::bayes($shortName:p) $parents
     set ::bayes($shortName:c) $children
     set ::bayes($shortName:w) $weights
     lappend ::bayes(nodes) $shortName
 }

 ###################################################################### 
 # start of main script below
 
 # do the breast cancer example
 addNode bc {} {s} 0.01 "Have Breast Cancer"
 addNode s {bc} {} {0.80 0.096} "Screened Positive"
 initBeliefs
 puts "Percent who have breast cancer given a positive screening: [infer +bc +s]"
 
 # now do the Tcl/Python example
 
 # first clear away the old Bayesian network
 array unset ::bayes
 # then populate it with nodes
 addNode kt {} {ll osg} 0.20 "Knows Tcl"
 addNode kp {} {ll osg} 0.25 "Knows Python"
 addNode ll {kt kp} {} {0.80 0.75 0.50 0.20} "Language Lawyer"
 addNode osg {kt kp} {h} {0.60 0.58 0.30 0.05} "OS Geek"
 addNode h {osg} {} {0.95 0.33} "Hireable"
 # initialize the beliefs network
 initBeliefs
 # now make queries
 puts "Percent of all programmers who are hireable: [evidence +h]"
 puts "Percent hireable who know Tcl and Python: [support +h +kt+kp]"
 puts "Percent of non-OS geeks who know just Python and not Tcl: [support ~osg +kp~kt]"
 puts "Percent who know Tcl but are not language lawyers: [infer +kt ~ll]"