#!/bin/sh
# the next line restarts using wish \
exec wish "$0" "$@"

# mutual.tcl - evaluates mutual info between x and y, you get to tweak px,Q
#                                   David J C MacKay (1999)

#
# set these two flags to 1 to get a stand-alone version
# as used by me; leave them at 0 to get plug-in netscape
# compatible version:
#
global verbose ; set verbose 0
global ownwindow ; set ownwindow 0
global version ; set version 2.1

set ownwindow 0
set verbose   0
if {$verbose>=1} {
    puts "Mutual Information $version              David J C MacKay          1999"
    puts "                                    written under tcl 8.0 on linux"
}
# changes
#########################################################################
# Feb 2000, added nice canvas
######################################################################
# contents:
#
######################################################################
global showenergytext ; set showenergytext 0
global showplaintext ; set showplaintext 1
global font ; set font "Courier 20 bold"

frame .th
pack  .th
set w .th
set top .
if {$ownwindow>=1} {
    wm title . "Mutual Info Version $version"
    wm iconname . "mutual"
    wm geometry . +10+10
}

#######################################################################
#
#  variables
#
#######################################################################
proc setTypes { } {
    global I J number atypes htypes ptypes color palecolor Active qtypes nicename probname postname
    # foreach input we need a type of probability
    
# set types "px q py"
# atypes = the types that can be tweaked with an a bar
# ptypes = the types that are functions
# htypes = probabilities that should have their h computed
    set atypes "px "
    set qtypes ""
    for {set i 0} {$i < $I} {incr i} {
	set qtypes "$qtypes q$i"
	set atypes "$atypes q$i"
	set number(q$i) J
    }
    set htypes "px $qtypes py"
    set ptypes "py"

    set color(px) "red"
    
    set color(q0) "springgreen1"
    set color(q1) "green2"
    set color(q2) "darkolivegreen2"
    set color(q3) "seagreen2"
    set color(q4) "green4"
    set palecolor(q0) "springgreen1"
    set palecolor(q1) "palegreen1"
    set palecolor(q2) "darkolivegreen2"
    set palecolor(q3) "seagreen2"
    set palecolor(q4) "palegreen4"
    set color(py) "blue"

    set palecolor(px) "pink"
    set palecolor(py) "skyblue"

    set palecolor(HXY)  "yellow1"
    set palecolor(HXgY) "pink2"
    set palecolor(HYgX) "skyblue2"
    set palecolor(IXY)  "gold1"

    set nicename(HXY)  "H(X,Y)"
    set nicename(HXgY) "H(X|Y)"
    set nicename(HYgX) "H(Y|X)"
    set nicename(IXY)  "I(X;Y)"
    set nicename(px)  "H(X)"
    set nicename(py)  "H(Y)"
    set nicename(q0) "H(Y|x=0)"
    set nicename(q1) "H(Y|x=1)"
    set nicename(q2) "H(Y|x=2)"
    set nicename(q3) "H(Y|x=3)"
    set nicename(q4) "H(Y|x=4)"

    set probname(px)  "P(x="
    set probname(py)  "P(y="
    set probname(q0) "Q(y="
    set probname(q1) "Q(y="
    set probname(q2) "Q(y="
    set probname(q3) "Q(y="
    set probname(q4) "Q(y="

    set postname(px)  ")"
    set postname(py)  ")"
    set postname(q0) "|0)"
    set postname(q1) "|1)"
    set postname(q2) "|2)"
    set postname(q3) "|3)"
    set postname(q4) "|4)"

    set number(px) I
    set number(py) J

    foreach t [concat $ptypes $atypes] {
	set Active($t) 1
    }
}

#######################################################################
#
#                      procedures
#
#######################################################################

global thewidth ;    set thewidth 50 ;# widths of sliders
global minwidth ;  set minwidth 8
global maxwidth ;  set maxwidth 50

proc setps { l t } {
    global ps p
    if {$p($l,$t)<1.0} { 
	set ps($l,$t) [string range $p($l,$t) 1 end]
    } else {
	set ps($l,$t) $p($l,$t) 
    }
}

proc makeEnergylevels { w energylevels } {
    global I Emin Emax ocmax ocmin betamin Z  T  J
    global energy types color meanoc verbose palecolor
    global thewidth  number atypes ptypes  htypes p minwidth maxwidth ps
    global showenergytext showplaintext probname postname

    # make our own local frame
    set e $energylevels.e
    catch { destroy $e }
    pack [frame $e]

    # redefine the slider widths
    set maybethewidth [expr {int(400.0/(($I+1)*($J+1)-1))}]
#    puts $maybethewidth
    if {$maybethewidth<$minwidth} {set maybethewidth $minwidth}
    if {$maybethewidth>$maxwidth} {set maybethewidth $maxwidth}
    set thewidth $maybethewidth 

    # make all the a bars
    foreach  t [concat $atypes] {
	set L [set $number($t)]
#	puts $t
#	puts $L
	for { set l 0 } { $l < $L } { incr l } {
	    set energy($l,$t) -4
	    catch { destroy $e.ee$l$t  ; destroy $e.e$l$t }
	    if {$verbose>=2} { puts "elmax = $Elmax" }
	    set ee$l [frame $e.ee$l$t] ; pack [set ee$l] -in $e -side left -pady 2
	    set e$l [scale $e.e$l$t -orient vertical -length 280  \
		    -from $Emin -to $Emax   -width $thewidth -sliderlength 8 -background $color($t) \
		    -borderwidth 0 -showvalue 0  \
		    -variable energy($l,$t) -tickinterval 0 -resolution 0.02  \
		    -bigincrement 0     -command {computeMicrostateProbs}]
	    pack [set e$l] -in [set ee$l] -side top -pady 2 -padx 2 -expand 1 -fill x

	    catch { destroy $e.et$l$t  }
	    if {$showenergytext} {
		set p$l$t [label $e.et$l$t -background $palecolor($t) \
			-textvariable energy($l,$t) -width 2 -anchor w]
		pack [set p$l$t] -in [set ee$l] -side top -pady 2   -padx 2 -expand 1 -fill x
	    }
	    
	    if {$showplaintext} {
		set p$l$t [label $e.et$l$t -background $palecolor($t) \
			-text "$probname($t)$l$postname($t)" -width 2 ]
		pack [set p$l$t] -in [set ee$l] -side top -pady 2   -padx 2 -expand 1 -fill x
	    }

	    # show probs too.
	    catch { destroy $e.p$l$t  }
	    set p$l$t [scale $e.p$l$t -orient vertical -length 80  \
		    -from $ocmax -to $ocmin   -width $thewidth -sliderlength 6 \
		    -borderwidth 0 -showvalue 0  -background \
		    $color($t) \
		    -variable p($l,$t) -tickinterval 0 \
		    -resolution 0.0001 ]
	    pack [set p$l$t] -in [set ee$l] -side top  -pady 2  -padx 2 -expand 1 -fill x
	    catch { destroy $e.pt$l$t  }
	    # feb 2000: replaced p(l,t) by ps(l,t) for "string"
	    set p$l$t [label $e.pt$l$t -background $palecolor($t) \
		    -textvariable ps($l,$t) -width 2 -anchor w]
	    pack [set p$l$t] -in [set ee$l] -side top  -pady 2  -padx 2 -expand 1 -fill x

	    catch { destroy $e.lpt$l$t  }
#	    set p$l$t [label $e.lpt$l$t -background $palecolor($t) \
#		    -textvariable logprob($l,$t) -width 2 -anchor w]
#	    pack [set p$l$t] -in [set ee$l] -side top  -padx 2  -pady 2 -expand 1 -fill x
	}
    }
#    raise    $elabels

    ##################################### 
    # insert the canvas here. put a "probability and canvas frame"
    set eepandc [frame $e.eepandc] ; pack [set eepandc] -in $e -side left -pady 2
    set eec     [frame $e.eec    ] ; pack [set eec    ] -in $eepandc -side top -pady 2

    global c black
    set black "#004444"
    set width 100
    set height 250
    set c $e.c
    catch {destroy $c}
    canvas $c -relief sunken -borderwidth 2 -width $width -height $height -background $black
#    set bg [lindex [$c config -bg] 4]
#    puts $c  ;
    pack $c -in $eec -side top 
    makecanvas $c 
    set eep     [frame $e.eep    ] ; pack [set eep    ] -in $eepandc -side bottom -pady 2
    #####################################
    # make all the p bars p(y)
    foreach  t [concat $ptypes] {
	set L [set $number($t)]
	for { set l 0 } { $l < $L } { incr l } {
	    set ee$l [frame $e.ee$l$t] ; pack [set ee$l] -in $eep -side left -pady 2
	    if {$showplaintext} {
		set p$l$t [label $e.et$l$t -background $palecolor($t) \
			-text "$probname($t)$l$postname($t)" -width 2 ]
		pack [set p$l$t] -in [set ee$l] -side top -pady 2   -padx 2 -expand 1 -fill x
	    }
	    # show probs too.
	    catch { destroy $e.p$l$t  }
	    set p$l$t [scale $e.p$l$t -orient vertical -length 80  \
		    -from $ocmax -to $ocmin   -width $thewidth -sliderlength 6 \
		    -borderwidth 0 -showvalue 0  -background \
		    $color($t) \
		    -variable ps($l,$t) -tickinterval 0 \
		    -resolution 0.0001 ]
	    pack [set p$l$t] -in [set ee$l] -side top  -pady 2  -padx 2 -expand 1 -fill x
	    catch { destroy $e.pt$l$t  }
	    set p$l$t [label $e.pt$l$t -background $palecolor($t) \
		    -textvariable p($l,$t) -width 2 -anchor w]
	    pack [set p$l$t] -in [set ee$l] -side top  -pady 2  -padx 2 -expand 1 -fill x

	}
    }


    catch { destroy $e.hframe }
#  entropies (these don't quite belong here... 
    set allzframe [frame $e.hframe]
    pack $allzframe -in $e -side right -padx 2
    pack [label $allzframe.zl -text "entropies"]  \
	    -side top -pady 0 -padx 10
    global H IXY
    global nicename
    foreach t [concat $htypes] {
	set H($t) 0.0 ;
	set zframe [frame $allzframe.z$t]
	set z$t [label $zframe.z  -width 6 -anchor w -background \
			$palecolor($t) -borderwidth 2 -textvariable H($t)]
	set zl$t [label $zframe.zl  -width 10 -anchor w -background \
			$palecolor($t) -borderwidth 2 -text "$nicename($t):"]
	pack [set zl$t] [set z$t]  -side left -pady 0 -padx 0
	pack [set zframe]  -side top -pady 6 -padx 6
    }

    foreach t [concat "HXY HXgY HYgX IXY"] {
	set zframe [frame $allzframe.z$t]
	set z$t [label $zframe.z  -width 6 -anchor w -background \
			$palecolor($t) -borderwidth 2 -textvariable $t]
	set zl$t [label $zframe.zl  -width 10 -anchor w -background \
			$palecolor($t) -borderwidth 2 -text "$nicename($t):"]
	pack [set zl$t] [set z$t]  -side left -pady 0 -padx 0
	pack [set zframe]  -side top -pady 6 -padx 6
    }


        # slider for IXY
    catch { destroy $e.iframe }
    set allzframe [frame $e.iframe]
    pack $allzframe -in $e -side right -padx 2
    foreach t [concat "IXY"] {
	set zframe [frame $allzframe.slider$t]
	set z$t  [scale $e.p$l$t -orient vertical -length 400  \
		    -from 1.2 -to 0.0   -width $maxwidth -sliderlength 10 \
		    -borderwidth 0 -showvalue 0  -background \
		    $palecolor($t) \
		    -variable $t -tickinterval 0.1 \
		    -resolution 0.00001 ]
	pack [set z$t]   -pady 0 -padx 0
	pack [set zframe]  -side top -pady 6 -padx 6
    }

}




# make a frame called l within frame controls, and associate these buttons 
# with an integer called L
proc adjustableInteger { w controls l Lname Lstring } {
    frame $w.$l 
    pack $w.$l -in $controls -side left  -pady 2  -padx 2 -anchor w
    button $w.$l.l -text "$Lstring:" -padx 0 -pady 0 -borderwidth 1 -command {}
    button $w.$l.up -text ">" -padx 0 -pady 0 -borderwidth 1 -command "incr $Lname"
    button $w.$l.dn -text "<" -padx 0 -pady 0 -borderwidth 1 -command "incr $Lname -1"
    bind  $w.$l.up  <3> "$w.$l.dn invoke"
    bind  $w.$l.dn  <3> "$w.$l.up invoke"
    entry $w.$l.n -textvariable $Lname -width 3 -borderwidth 1
    pack $w.$l.l $w.$l.dn $w.$l.up $w.$l.n -in $w.$l -side left

}

proc makeControls { w controls } {
    pack [button $w.restart -text "Restart" -borderwidth 1 \
	    -command {restart}] -in $controls  -side left \
	    -pady 0 -padx 2 -anchor w
    adjustableInteger $w $controls "i" "II" "I"
    adjustableInteger $w $controls "j" "JJ" "J"
    adjustableInteger $w $controls "ema" "Emin" "Emin"
    adjustableInteger $w $controls "emi" "Emax" "Emax"
#    adjustableInteger $w $controls "w" "thewidth" "width"
}

#
# packing procedures
#
proc restart { } {
    global N L II I JJ J
    upvar w w   energylevels energylevels  controls controls  
    upvar microstates microstates

    set I $II
    set J $JJ
    global ocmin ; set ocmin 0 ;
    global ocmax ; set ocmax 1 ;

    setTypes
    makeEnergylevels $w $energylevels
}

# invoked when the energy level scales are touched
# finds probabilities 
proc computeMicrostateProbs { {junk 0} } {
    global energy I J T atypes  p  Z KL KLqp logprob ps
    global Active number htypes
    foreach  t [concat $atypes] {
	set L [set $number($t)]
	set Z($t) 0.0
	for { set l 0 } { $l < $L } { incr l } {
	    if {$energy($l,$t)<0.0} {
		set pn($l,$t) [expr  exp(-$energy($l,$t)*log(2.0))]
	    } else {
		set pn($l,$t) 0.0
	    }
	    set Z($t) [expr $Z($t) + $pn($l,$t)]
	}
	for { set l 0 } { $l < $L } { incr l } {
	    set p($l,$t) [expr   $pn($l,$t)/$Z($t)]
#	    set logprob($l,$t) [expr   -log($p($l,$t))/log(2.0)]
	    setps $l $t
	}
    }

    computeDerivedProbs
    
    global H htypes
    # compute entropies
    foreach  t [concat $htypes] {
	set L [set $number($t)]
	set H($t) 0.0
	for { set l 0 } { $l < $L } { incr l } {
	    if {$p($l,$t)>0.0} {
		set H($t) [expr $H($t) - $p($l,$t) * log($p($l,$t))]
	    }
	}
	set H($t) [expr $H($t)/log(2.0)]
    }
    
    global HXY IXY HXgY HYgX qtypes p I
    set HYgX 0.0
    for {set i 0} {$i < $I} {incr i} {
	set HYgX [expr ($HYgX +	$p($i,px) * $H(q$i))]
    }
    set HXY  [expr $H(px) + $HYgX]
    set HXgY [expr $HXY - $H(py)]
    set IXY  [expr $H(px) - $HXgY]
}
# invoked when the energy level scales are touched
# finds probabilities 
proc computeDerivedProbs { {junk 0} } {
    global energy I J T ptypes  p  Z KL KLqp logprob
    global Active number
    foreach  t [concat $ptypes] {
	set L [set $number($t)]
	# L = number of outputs, M number of inputs
	set M [set $number(px)]
	set Z($t) 0.0
	for { set l 0 } { $l < $L } { incr l } {
	    set sum 0
	    for { set m 0 } { $m < $M } { incr m } {
		set sum [expr ($sum + $p($m,px) * $p($l,q$m))]
	    }
	    set p($l,$t) [expr $sum]
	    setps $l $t
	}
    }
##################### new in feb 2000: canvas stuff
    global c wfactor color black
    foreach  t [concat $ptypes] {
	set L [set $number($t)]
	# L = number of outputs, M number of inputs
	set M [set $number(px)]
	for { set l 0 } { $l < $L } { incr l } {
	    for { set m 0 } { $m < $M } { incr m } {
		set width [expr ($wfactor * $p($l,q$m))]
		if {$p($l,q$m)<0.0001} {
		    $c itemconfig edge$m$l -fill $black
		} else {
		    $c itemconfig edge$m$l -fill [set color(q$m)] 
		}
		$c itemconfig edge$m$l -width $width
	    }
	}
    }
    
}
global wfactor ; set wfactor 10 ; 

####################################################################
#
# set up windows 
# 
####################################################################

global II ; set II 2
global JJ ; set JJ 2
# T = number of transitions
global Emin ; set Emin -8 ;
global Emax ; set Emax 1 ;# highest energy for a level

bind . <Control-x> "destroy ."
bind . <Control-q> "destroy ."
bind . <q> "destroy ."
bind . <Control-c> "destroy ."
bind . <Control-z> "destroy ."


global c 
# c is the canvas

proc makecanvas {c} {
    global I J number atypes htypes ptypes color palecolor Active qtypes nicename probname postname
    global font
    set leftx 30
    set rightx 70
    set leftxl 18 ;# label locations
    set rightxl 85
    set dy   30
    set topy 10
#    catch {$c delete text}
#    catch {$c delete edge}

    for {set i 0 } {$i < $I} {incr i} {
	set yi [expr {$dy * $i + $topy}]
	for {set j 0} {$j < $J} {incr j} {
	    set yj [expr {$dy * $j + $topy}]
	    set item [$c create line ${leftx} ${yi} ${rightx} ${yj} \
		-width 2 -fill [set color(q$i)] -tags [list edge edge$i$j ]]
	}
    }
    for {set i 0 } {$i < $I} {incr i} {
	set yi [expr {$dy * $i + $topy}]
	set item [$c create text ${leftxl}  ${yi} -anchor e \
		-text "$i" -font $font -tags text -fill indianred1]
    }
    for {set i 0 } {$i < $J} {incr i} {
	set yi [expr {$dy * $i + $topy}]
	set item [$c create text ${rightxl}  ${yi} -anchor w \
		-text "$i" -font $font -tags label -fill skyblue1]
    }
}



frame $w.controls
set controls $w.controls
pack $controls -side top
makeControls $w $controls

frame $w.middlerow
set middlerow $w.middlerow
pack $middlerow -side top

frame $w.energylevels
set energylevels $w.energylevels
pack $energylevels -in $middlerow -side top -expand y -fill x

# buttons to do with overall control (quit, etc.)
frame $w.buttons
pack $w.buttons -side bottom -fill x -pady 2 -expand 1

#
# overall control buttons
#
button $w.dismiss -text Quit -command "destroy $top"
button $w.help -text Help -command "help"

pack $w.dismiss  $w.help \
	-in $w.buttons -side left -fill x -expand 1 -anchor w  -padx 3 -pady 1
############
#  end bottom row
############

# make it happen!
restart

#####################################################################

proc help { } {
    set w .help
    catch {destroy $w}
    global ownwindow
    if {$ownwindow>=1} {
	toplevel $w
	wm geometry $w +10+10
	bind $w <Control-c> "destroy $w"
    } else {
	pack [frame $w]
	bind $w <Control-c> "destroy $w"
    }
    frame $w.buttons
    pack $w.buttons -side bottom -fill x -pady 2 -expand 1 -padx 4
    button $w.buttons.dismiss -text Dismiss -command "destroy $w"
    pack $w.buttons.dismiss  -side left -fill x -expand 1 -padx 4

  text $w.t -background white -height 24 -wrap word\
            -xscrollcommand "$w.xscroll set" \
            -yscrollcommand "$w.yscroll set" \
            -setgrid 1 -highlightthickness 0 -pady 2 -padx 3
        scrollbar $w.xscroll -command "$w.t xview" \
            -highlightthickness 0 -orient horizontal
        scrollbar $w.yscroll -command "$w.t yview" \
            -highlightthickness 0 -orient vertical

pack $w.yscroll -side right -fill y
pack $w.xscroll -side bottom -fill x
pack $w.t -expand yes -fill both

    $w.t insert 0.0 \
{Mutual Information        - author David J C MacKay  mackay@mrao.cam.ac.uk

 Make your own channel, and find its mutual information.
 Things to notice:
      I(X;Y) is a convex function of Q(y|x)
      I(X;Y) is a concave function of P(x)
     
 (I(X;Y) may appear to fluctuate when close to zero, because scientific notation kicks in)

 Adjust the adjustable probabilities by yanking the upper red and green sliders.

 General layout:
     Red stuff: input distribution P(x) 
     Green stuff: conditional distribution Q(y|x)
     Blue stuff:  output distribution P(y) = sum_x P(x) Q(y|x)

 Shortcuts: 
     C-r     reset

 Maximum number of inputs and outputs is 5. Recommended settings are 2,3,4.
}
}