#!/usr/local/bin/wish -f
#
#  LearnHMM - a program for learning HMM
#
#  coded by Shigeki Sagayama, JAIST, 26 Oct 1998

#------------------------------------------------------------------------------
# help commands
#------------------------------------------------------------------------------
proc helpIntro {} { helpDisp {help: 初めに

          北陸先端科学技術大学院大学 情報科学研究科   嵯峨山茂樹
             sagayama@jaist.ac.jp, http://www-ks.jaist.ac.jp/

このプログラムはＨＭＭ(hidden Markov model, 隠れマルコフモデル)の学習を助ける
ためのソフトウェアです。目的は、

(1) 確率モデルとしてＨＭＭの定義を理解する。モデルパラメータを自由に変えられる
    ようにしてあります。

(2) 非定常な時系列を生成するモデルとしてＨＭＭを理解する。このために、確率的に
    信号系列を生成し、その確率を表示するようにしてあります。

(3) 与えられた信号に対する確率評価機能としてＨＭＭを理解する。このために、任意
    の信号系列を入力できるようにしてあります。

プログラム言語としてはグラフィカルなスクリプト言語である Tcl/Tk を用いています。
Tcl/Tk は GUI を活用したインタラクティブなシステムを作るのに適しているのですが、
ＨＭＭのような確率計算を記述するには向いていないので、プログラミングには苦労し
ます。しかし、無料のソフトで、広く普及していて、インタラクティブなシステムを作
るのに適している言語は Tcl/Tk 以外にないので、処理速度が遅い点は我慢して貰うこ
とにして使っています。

以上。}}
#------------------------------------------------------------------------------
proc helpCopyright {} { helpDisp {\
Help: Copyright

◆ このプログラムの著作権は、作者・嵯峨山茂樹 (sagayama@jaist.ac.jp) に属し
ます。

作者は、教材作りを目的に Tcl/Tk を勉強し始めてまだ一ヶ月ほどなので、プログラム
の書き方が下手な箇所は随所に見られることと思います。また、ＨＭＭのようなパター
ン認識計算分野を Tcl/Tk のスクリプトで行なうという発想がそもそも特異なので、改
善の余地は多くあることと思います。改善の御意見を歓迎します。

不完全な代物ではありますが、しかし、作者に無断で改善・改造して、他の人に配布す
ることは避けて下さい。そこで、以下に使用条件を述べます。

◆ このソフトウェアを改変することなくそのまま全部を再配布することは許可しま
す。大いに使って頂いて結構です。ただし、改善・改変して再配布しないで下さい。
この著作権に関する宣言の部分も含め、一字一句変えずにこのままならば再配布は自
由です。但し、漢字コードの変換などは構いません。

◆ 作者の許可なく改変(一部分のみの使用を含む)したものは、再配布することを禁
じます。

◆ バグやプログラミングに関する御意見は、作者に御連絡下さい。バグ情報、改善意
見、改良の努力、利用経験、応用情報などは歓迎致します。

以上。}
}
#------------------------------------------------------------------------------
proc helpModel {} { helpDisp {help: 確率モデルの操作

モデルのパラメータ操作は次の２つのモードで行なえます。
  (1) 状態入力   .......   ＨＭＭ状態出力確率分布の平均・分散を変えるモード
  (2) 遷移入力   .......   ＨＭＭ状態間の遷移確率を変えるモード

３種の入力モードの切替え操作は、つぎのいずれかを用いて行ないます。
  (1) 入力モード表示ボタンが目的のモードになるまでクリックを繰り返す。
  (2) ショートカットキー操作
        d  ...  信号系列入力
        s  ...  状態出力分布入力
        t  ...  遷移確率入力

状態入力モードでの操作
  (1) マウス左ボタンをクリック・ドラッグ
      分布の分散値を変化させます。
  (2) マウス中ボタンをクリック・ドラッグ
      分布の平均ベクトルを変化させます。

遷移入力モードでの操作
  (1) マウス左ボタンで矢印をクリック
      どの状態間遷移の確率を変えるかを指定します。
  (2) マウス左ボタンで「遷移」スケールをドラッグ
      指定された状態間遷移の確率を変えます。

以上。}}
#------------------------------------------------------------------------------
proc helpInput {} { helpDisp {help: 信号系列の入力

信号系列の入力は次の２種があります。
  (1) 手動データ入力   .....   データ系列を画面から与えるモード
  (2) 自動生成入力   .......   データを自動的にＨＭＭから生成

３種の入力モードの切替え操作は、つぎのいずれかを用いて行ないます。
  (1) 入力モード表示ボタンが目的のモードになるまでクリックを繰り返す。
  (2) ショートカットキー操作
        d  ...  データ系列入力
        s  ...  状態出力分布入力
        t  ...  遷移確率入力

信号入力モードでの操作
  (1) マウス左ボタンをクリック
      クリックした場所に既にデータが存在する場合は、「選択」を意味します。
      既にデータが存在しない場合は、新たにデータを追加します。
  (2) マウス中ボタンをクリック・ドラッグ
      データを移動します。
  (3) マウス右ボタンをクリック
      該当するデータを消去します。
  (4) キーボード <c> をタイプ
      「選択」されているデータをコピーして少しずらしたデータを系列に挿
      入します。
  (5) 「信号長」ボタンをクリック
      信号をクリアします。長さは０になります。
  (6) 「生成」をクリック
      表示されているＨＭＭに基づいて、信号系列をランダムに生成します。生成され
      る確率密度が log prob = .... に表示され、同時に右下に赤い棒の長さで表さ
      れます。

以上。}}
#------------------------------------------------------------------------------
proc helpViterbi {} { helpDisp {help: Viterbi 経路探索

与えられた信号系列に対し、ＨＭＭを最適にマッチングし、その状態遷移経路を求め、
表示します。

  (1) log prob = .... と表示されているボタンをクリック。
  (2) 信号系列に状態との対応づけ(Viterbi alignment)がなされて、色で表示されま
      す。
  (3) ボタンに、log prob = .... と Viterbi スコアが表示されます。同時に右下に
      赤い棒の長さで表されます。

以上。}}
#------------------------------------------------------------------------------
proc helpModels {} { helpDisp {help: モデルの切替え

モデルを退避し、また復元することができます。

  (1) ボタン「モデル」をクリックするとポップアップメニューが表示されます。
  (2) 「モデルを記憶」をクリックすると、現在のモデルが番号をつけて記憶されます。
  (3) 以後、「モデル(n)復元」をクリックすると、n 番目に記憶したモデルが復元さ
      れて表示され、計算に用いられます。

ひとつのモデルは、音声認識では単語音声のモデルに相当します。同一の入力系列に対
して、複数のモデルにより Viterbi 経路と尤度を計算し、比較して、最も尤度高いモ
デルに対応する単語を認識結果にするのが、音声認識の原理です。

以上。}}
#------------------------------------------------------------------------------
proc helpLearn {} { helpDisp {help: 学習のヒント

このソフトウェアを用いてＨＭＭの学習をしましょう。たとえば、つぎのよう
なＨＭＭの性質を実際に試してみるとよいでしょう。

(1) 何度も信号生成を行なってみる。信号系列長、状態継続時間、サンプルベクトル値
    がいかにばらつくか。

(2) 遷移確率を変えて、信号生成を行なって見る。状態継続時間がどのように変化する
    か。

(3) 状態出力確率分布の平均ベクトルを変えてみる。出力信号の分布がどのように変化
    するか。

(4) 状態出力確率分布の分散を変えてみる。出力信号の分布がどのように変化するか。

(5) 任意の信号系列を入力して、Viterbi 尤度を計算する。どのように時間の整合が取
    られるか。信号系列とモデルの状態系列がかなり離れている場合に、Viterbi 経路
    と尤度はどのようになるか。

(6) 単語音声認識の原理を理解する。同一の信号系列に対して、複数のモデルによる 
    Viterbi 尤度を計算し、比較して、最も尤もらしいモデルを選ぶ。これが認識の原
    理である。

以上。}}
#------------------------------------------------------------------------------
### Create a new message window to show the specified contents ###
proc helpDisp { contents } {
  set w .basic
  catch {destroy $w}
  toplevel $w
  wm title $w "Fourier Transform"
  wm iconname $w "Fourier Transform"
  button $w.ok -text OK -command "destroy $w" -fg black -bg yellow
  text $w.t -relief raised -bd 2 -yscrollcommand "$w.s set" -setgrid true \
    -fg black -bg white
  scrollbar $w.s -relief flat -command "$w.t yview"
  pack $w.ok -side bottom -fill x
  pack $w.s -side right -fill y
  pack $w.t -expand yes -fill both
  $w.t insert 0.0 $contents
  $w.t mark set insert 0.0
  bind $w <Any-Enter> "focus $w.t"
}
#------------------------------------------------------------------------------
# initial values
#------------------------------------------------------------------------------
set W 801
set H 600
set H0 150
set Sy 100

# initial number of states ( arbitrary )
set Ns 3

# changeable
set PrevN 0
set Np 0
set Nm 0
set PointList { }
set This 1
set Smooth off
set cbg #ffe1bc

set numb 0; set total 0; set seconds 0; set hundredths 0;
set Delta 100

# initialize $Mode = "data",  (others: "state", "trans")
set Mode "state"

### colors and initial HMM parameters ###
proc initHMM {} {
  global W H H0 Ns Fg Bg Cx Cy Vx Vy Trans
  for { set i 1 } { $i <= $Ns } { incr i } {
    set fg [format "#%02x%02x%02x" \
           [expr 156+int(100*cos(3.14159*2*(double($i)/$Ns)))] \
           [expr 156+int(100*cos(3.14159*2*(double($i)/$Ns+1.0/3)))] \
           [expr 156+int(100*cos(3.14159*2*(double($i)/$Ns+2.0/3)))] ]
    set bg [format "#%02x%02x%02x" \
           [expr 128+int(48*cos(3.14159*2*(double($i)/$Ns)))] \
           [expr 128+int(48*cos(3.14159*2*(double($i)/$Ns+1.0/3)))] \
           [expr 128+int(48*cos(3.14159*2*(double($i)/$Ns+2.0/3)))] ]
    set j state$i
    set Fg($j) $fg
    set Bg($j) $bg
    set Cx($j) [expr $i*$W/($Ns+1)+$W*((2*$i)%$Ns/($Ns*10.))]
    set Cy($j) [expr $H0+(2111*$i)%($H-$H0)]
    set Vx($j) [expr $W*($i*247)%63+50]
    set Vy($j) [expr $H*($i*391)%43+30]
    set Trans($j) [expr ($i*17)%29/100.+0.3]
  }
}
initHMM
#------------------------------------------------------------------------------
# display setup
#------------------------------------------------------------------------------
### window title ###
wm title . "ShapeDesign"

### message ###
message .m -justify left -relief raised -bd 2 -width [expr $W+30] \
  -bg navy -fg cyan -text \
"LearnHMM  (sagayama@jaist.ac.jp) - 試作中"

#### canvas ####
canvas .c -width $W -height $H -bg wheat
#------------------------------------------------------------------------------
### buttons: start, stop, quit, time, numb ###
frame .f

button .f.generate -text "生成" -command {
  genPoints
} -bg slateblue -fg white

button .f.numb -text "信号長 0" -command {
  erasePoints
} -width 10 -bg seagreen -fg white
proc showNp { } {
  global Np
  .f.numb config -text [format "信号長 %d" $Np]
}

# button .f.curve -text "曲線／直線" -command {
#   if { $Smooth } { set Smooth 0 } else { set Smooth 1 }
#   drawShape
# } -bg orangered -fg white

button .f.mode -text "状態モード(s)" -command {
  switch $Mode {
    "data" { set Mode "state" }
    "state" { set Mode "trans" }
    "trans" { set Mode "data" }
  }
  showMode
} -bg maroon -fg white
proc showMode {} {
  global Mode
  switch $Mode {
    "data"  { .f.mode configure -text "信号モード(d)" }
    "state" { .f.mode configure -text "状態モード(s)" }
    "trans" { .f.mode configure -text "遷移モード(t)" }
  }
}

button .f.stop -text "出力" -command {
# simple output: (OK)
  puts stdout [.c coord $Shape]
} -bg slateblue -fg white

menubutton .f.nst -text "状態数" -menu .f.nst.m -relief raised \
  -fg #ffffff -bg #cd2800
menu .f.nst.m -fg #ffffff -bg #910d2b
.f.nst.m add command -label "1 状態" -command { changeNs 1 }
.f.nst.m add command -label "2 状態" -command { changeNs 2 }
.f.nst.m add command -label "3 状態" -command { changeNs 3 }
.f.nst.m add command -label "4 状態" -command { changeNs 4 }
.f.nst.m add command -label "5 状態" -command { changeNs 5 }
.f.nst.m add command -label "6 状態" -command { changeNs 6 }
.f.nst.m add command -label "7 状態" -command { changeNs 7 }
.f.nst.m add command -label "8 状態" -command { changeNs 8 }
.f.nst.m add command -label "9 状態" -command { changeNs 9 }
.f.nst.m add command -label "10 状態" -command { changeNs 10 }
.f.nst.m add command -label "11 状態" -command { changeNs 11 }
.f.nst.m add command -label "12 状態" -command { changeNs 12 }
.f.nst.m add command -label "13 状態" -command { changeNs 13 }
.f.nst.m add command -label "14 状態" -command { changeNs 14 }
.f.nst.m add command -label "15 状態" -command { changeNs 15 }
.f.nst.m add command -label "16 状態" -command { changeNs 16 }
.f.nst.m add command -label "17 状態" -command { changeNs 17 }
.f.nst.m add command -label "18 状態" -command { changeNs 18 }
.f.nst.m add command -label "19 状態" -command { changeNs 19 }
.f.nst.m add command -label "20 状態" -command { changeNs 20 }
### restart this program ###
proc changeNs { ns } {
  global Ns Nm
  erasePoints
  .c delete all
  set Ns $ns
  set Nm 0
  initHMM
  drawHMM
  drawInitStates
  drawGrid 100
  initbarScore
  update idletasks
}

### model stack ###
menubutton .f.model -text "モデル" -menu .f.model.m -relief raised \
  -fg #ffffff -bg #a37100
menu .f.model.m -fg #ffffff -bg #a8371c
.f.model.m add command -label "モデルを記憶" -command {
  saveModel
}
### save models in the memory ###
proc saveModel { } {
  global Nm Cx Cy Vx Vy MemCx MemCy MemVx MemVy
  incr Nm
  eval .f.model.m add command -label "モデル($Nm)復元" \
    -command \{ recallModel $Nm \}
  foreach dist [.c find withtag "state"] {
    set tags [lindex [.c itemconfigure $dist -tags] 4]
    set i [lindex $tags 0]
    set k $i/$Nm
#    puts $k
    set MemCx($k) $Cx($i)
    set MemCy($k) $Cy($i)
    set MemVx($k) $Vx($i)
    set MemVy($k) $Vy($i)
  }
}
### recall models from the memory ###
proc recallModel { id } {
  global Nm Cx Cy Vx Vy Sx Sy MemCx MemCy MemVx MemVy
#  puts "recall model $id"
  foreach dist [.c find withtag "state"] {
    set tags [lindex [.c itemconfigure $dist -tags] 4]
    set i [lindex $tags 0]
    set k $i/$id
    set Cx($i) $MemCx($k)
    set Cy($i) $MemCy($k)
    set Vx($i) $MemVx($k)
    set Vy($i) $MemVy($k)
    .c coords $dist \
      [expr $Cx($i)-$Vx($i)] [expr $Cy($i)-$Vy($i)] \
      [expr $Cx($i)+$Vx($i)] [expr $Cy($i)+$Vy($i)]
    set tag to$i
    .c coords [.c find withtag $tag] $Sx($i) $Sy $Cx($i) $Cy($i)
  }
}

button .f.score -text "log prob = (確率・尤度)" -relief raised \
  -bg #f5f037 -fg #000000 -command {
  Viterbi
}

menubutton .f.help -text "使用法" -menu .f.help.m -relief raised \
  -fg #ffffff -bg #32850d
menu .f.help.m -fg #ffffff -bg #1e6e1c
.f.help.m add command -label "初めに" -command helpIntro
.f.help.m add command -label "Copyright" -command helpCopyright
.f.help.m add command -label "モデル操作／変更" -command helpModel
.f.help.m add command -label "信号入力／生成" -command helpInput
.f.help.m add command -label "Viterbi経路探索" -command helpViterbi
.f.help.m add command -label "モデルの切替え" -command helpModels
.f.help.m add command -label "使い方のヒント" -command helpLearn

button .f.quit -text "終了" -command {exit} -bg #be0d2b -fg white
#------------------------------------------------------------------------------
### scales ###
frame .s

# Get the state transition probability
# scale .s.trans -label "遷移" -from 100 -to 0
scale .s.trans -from 100 -to 0 \
  -length [expr $H*2/3] \
  -orient vertical -command getTrans -bg grey30 -fg white
pack .s.trans -side top
proc getTrans value {
  global This Prob Trans Loop
  set tags [lindex [.c itemconfigure $This -tags] 4]
  if { [lsearch $tags "trans"] >=0 } {
    set i [lindex $tags 0]
    set Trans($i) [expr [.s.trans get]/100.]
    .c itemconfigure $Prob($i) -text $Trans($i)
    .c itemconfigure $Loop($i) -text [expr 1-$Trans($i)]
  }
}
.s.trans set 50

# Redraw the grids
# scale .s.delta -label "目盛" -from 100 -to 1
scale .s.delta -from 100 -to 1 \
  -length [expr $H/3] -bg darkgrey -fg white -orient vertical -command getDelta
proc getDelta value {
  eval .c delete [.c find withtag grid]
  drawGrid [.s.delta get]
}
.s.delta set 100
pack .s.delta -side top
#------------------------------------------------------------------------------
#  packing and binding
#------------------------------------------------------------------------------
pack .m -side top -fill both -expand yes
# pack .f.generate .f.numb .f.mode .f.stop .f.model .f.score \
#    .f.help .f.quit \
#   -side left -fill both -expand yes
pack .f.generate .f.numb .f.mode .f.nst .f.model .f.score \
   .f.help .f.quit \
  -side left -fill both -expand yes
pack .f -side top -fill both -expand yes
pack .s -side right
pack .c -side bottom -fill both -expand yes

# dangeous - remove the following two lines later
bind . <Control-c> {destroy .}
bind . <q> {destroy .}

focus .
#------------------------------------------------------------------------------
#  mouse actions
#------------------------------------------------------------------------------
bind .c <Button-1> {
  global Point Np This ShowN PointList
  set tags [lindex [.c itemconfigure current -tags] 4]
  switch $Mode {
    "data" {
      if { [lsearch $tags current] >=0 && [lsearch $tags "data"] >= 0 } {
      # existing point: focus on it
        .c itemconfigure current -fill red
         set curr [.c find withtag current]
        if { $This != $curr } { .c itemconfigure $This -fill blue }
        set This $curr
      } else {
      # no point existing here: generate one
        set x %x
        set y %y
        set This [addPoint $x $y red]
      }
    }
    "state" {
      set curr [.c find withtag current]
      set This $curr
    }
    "trans" {
      if { [lsearch $tags "trans"] >= 0 } {
        set This [.c find withtag current]
        .c itemconfigure trans -fill black
        .c itemconfigure current -fill red
        set i [lindex $tags 0]
        .s.trans config -bg $Bg($i)
        .s.trans set [expr int($Trans($i)*100)]
      }
    }
  }
}

# change the distribution variances
bind .c <B1-Motion> {
  global This Np ShowN
  if { $Mode == "state" } {
    set tags [lindex [.c itemconfigure current -tags] 4]
    if { [lsearch $tags "state"] >= 0 } {
      set i [lindex $tags 0]
      set x %x
      set y %y
      set Vx($i) [expr abs($x-$Cx($i))]
      set Vy($i) [expr abs($y-$Cy($i))]
      .c coords current \
         [expr $Cx($i)-$Vx($i)] [expr $Cy($i)-$Vy($i)] \
         [expr $Cx($i)+$Vx($i)] [expr $Cy($i)+$Vy($i)]
    }
  }
}

# move a point
bind .c <Button-2> {
  global This Np ShowN
  set tags [lindex [.c itemconfigure current -tags] 4]
#  puts stdout Button-2:$tags
  if { [lsearch $tags $Mode] >= 0 } {
    if { $Mode == "data" } { .c itemconfigure current -fill red }
    set curr [.c find withtag current]
    if { $This != $curr && $Mode == "data" } {
      if { [lsearch [.c itemconfigure $This -tags] "data"] >= 0 } {
        .c itemconfigure $This -fill blue
      }
    }
    set This $curr
  }
}

# drag a point
bind .c <B2-Motion> {
  global This Np ShowN
  if { [lsearch $tags $Mode] >= 0 } {
    if { $Mode == "data" } {
      set xx [expr [lindex [.c coords $This] 2]-5]
      set yy [expr [lindex [.c coords $This] 3]-5]
      set x %x
      set y %y
      .c move $This [expr $x-$xx] [expr $y-$yy]
      .c itemconfigure $This -tags [list $x $y data]
    } else {
      if { $Mode == "state" } {
        set tags [lindex [.c itemconfigure $This -tags] 4]
        set i [lindex $tags 0]
        set xx [expr [lindex [.c coords $This] 2]-$Vx($i)]
        set yy [expr [lindex [.c coords $This] 3]-$Vy($i)]
        set Cx($i) %x
        set Cy($i) %y
        .c move $This [expr $Cx($i)-$xx] [expr $Cy($i)-$yy]
        .c coords [.c find withtag to$i] $Sx($i) 100 $Cx($i) $Cy($i)
      }
    }
  }
}

# delete a point
bind .c <Button-3> {
  global Point Np This ShowN PointList
  set tags [lindex [.c itemconfigure current -tags] 4]
  if { [lsearch $tags data] >= 0 && [lsearch $tags current] >=0 } {
    set curr [.c find withtag current]
    delPoint $curr
    drawShape
  }
}

# release a button
bind .c <ButtonRelease-1> { drawShape }
bind .c <ButtonRelease-2> { drawShape }
bind .c <ButtonRelease-3> { }

# duplicate
bind . <c> { set This [dupPoint $This]; drawShape }

# delete
bind . <x> {
  if { $Mode == "data" } delPoint $This
  drawShape
}

# modify data
bind . <d> {
  set Mode "data"
  showMode
}
# modify state distributions
bind . <s> {
  set Mode "state"
  showMode
}
# modify transition probabilities
bind . <t> {
  set Mode "trans"
  showMode
}
#------------------------------------------------------------------------------
# grid
#------------------------------------------------------------------------------
### draw grid lines in grey ###
proc drawGrid { wid } {
  global W H H0
  for { set x 0 } { $x < $W } { incr x $wid } {
    set id [.c create line $x $H0 $x $H -fill darkgrey -tags grid]
    .c lower $id
  }
  for { set y $H0 } { $y < $H } { incr y $wid } {
    set id [.c create line 0 $y $W $y -fill darkgrey -tags grid]
    .c lower $id
  }
}
#------------------------------------------------------------------------------
# erase, add, delete, insert, duplicate, etc., data points
#------------------------------------------------------------------------------
### erase all data ###
proc erasePoints { } {
  global Np PointList
  eval .c delete [.c find withtag data]
  eval .c delete [.c find withtag line]
  set PointList {}
  set Np 0
  showNp
}
### add a point (click button 1 anywhere) ###
proc addPoint { x y color } {
  global Np PointList This
  if { [lsearch [lindex [.c itemconfigure $This -tags] 4] "data"] >= 0 } {
#    .c itemconfigure $This -fill $color
  }
  set id [.c create rectangle [expr $x-5] [expr $y-5] [expr $x+5] [expr $y+5]\
    -fill red -tags [list $x $y "data" "move"] ]
  lappend PointList $id
  incr Np
  set This $id
  .c itemconfigure $This -fill $color
  showNp
  return $id
}

### duplicate a point (click button 1 on an existing point)  ###
proc dupPoint { id } {
  global Np PointList This
  set tags [lindex [.c itemconfigure $id -tags] 4]
  set x [lindex $tags 0]
  set y [lindex $tags 1]
  set i [lsearch $PointList $id]
  set dup [.c create rectangle \
    [expr $x-10] [expr $y-10] [expr $x] [expr $y] \
    -fill red -tags [list [expr $x-5] [expr $y-5] "data"]]
  set PointList [linsert $PointList $i $dup]
  incr Np
  showNp
  return $dup
}

### delete a point (click button 3 on an existing point) ###
proc delPoint { id } {
  global Np PointList This
  set Np [llength $PointList]
  if { $Np <= 0 } { return 0 }
  for { set i 0 } { $i < $Np } { incr i } {
    if { [lindex $PointList $i] == $id } break
  }
  .c delete $id
  set PointList [concat [lrange $PointList 0 [expr $i-1]] \
                        [lrange $PointList [expr $i+1] end]]
  incr Np -1
  if { $id == $This } {
    if { $i >= $Np } { incr i -1 }
    if { $i < 0 } { set i 0 }
    set This [lindex $PointList $i]
    .c itemconfigure $This -fill red
  }
  if { $i < $Np } { return [lindex $PointList $i] } else { return 0 }
  showNp
}
### draw or redraw the time sequence shape ###
set Shape { }
proc drawShape { } {
  global Shape PointList PrevN Smooth Np
#  set Np [llength $PointList]
#  puts Np=$Np
  eval .c delete [.c find withtag line]
  if { $Np > 1 } {
    set pnts { }
    for { set i 0 } { $i < $Np } { incr i } {
      set point [lindex $PointList $i]
      set tags [lindex [.c itemconfigure [lindex $PointList $i] -tags] 4]
      set pnts [concat $pnts [lrange $tags 0 1]]
    }
#    puts $pnts
    if { [llength $pnts] > 0 } {
      set Shape [eval .c create line $pnts -fill black -smooth $Smooth \
        -tags line]
    }
  }
  set PrevN $Np
  showNp
}
#------------------------------------------------------------------------------
# draw an HMM topology and define transition probabilities
#------------------------------------------------------------------------------
proc drawHMM { } {
  global W H Ns Sx Sy Loop Prob Trans Fg Bg
  set y $Sy
  set dx [expr $W/($Ns+3)]
  .c create line [expr $dx+20] $y [expr $dx*2-20] $y -width 2 -tags a01
  .c create line [expr $dx*2-30] [expr $y-5] [expr $dx*2-20] $y \
    [expr $dx*2-30] [expr $y+5] -width 2 -tags arrow
  for { set i 1 } { $i <= $Ns } { incr i } {
    set x [expr $dx*($i+1)]
    set Sx(state$i) $x
    .c create oval [expr $x-20] [expr $y-20] [expr $x+20] [expr $y+20] \
      -fill $Bg(state$i) -width 2 -tags state$i
    .c create line [expr $x-10] [expr $y-17] [expr $x-20] [expr $y-50] \
      $x [expr $y-80] [expr $x+20] [expr $y-50] [expr $x+10] [expr $y-17] \
      -width 2 -smooth on -tags state$i
    .c create line [expr $x-18] [expr $y-26] [expr $x-10] [expr $y-17] \
      [expr $x-8] [expr $y-29] -width 2 -tags arrow
    eval set Loop(state$i) \
      [.c create text $x [expr $y-75] -text [expr 1-$Trans(state$i)] \
        -anchor s]
    eval set Prob(state$i) \
      [.c create text [expr $x+$dx/2] [expr $y-10] -text $Trans(state$i) \
      -anchor s]
    eval .c create line [expr $x+20] $y [expr $x+$dx-20] $y -width 2 \
      -tags \{ state$i trans \}
    .c create line [expr $x+$dx-30] [expr $y-5] [expr $x+$dx-20] $y \
      [expr $x+$dx-30] [expr $y+5] -width 2 -tags arrow
  }
}
drawHMM
#------------------------------------------------------------------------------
# draw output probability distributions associated with hidden states
#------------------------------------------------------------------------------
proc drawInitStates { } {
  global Ns Cx Cy Vx Vy Sx Sy Bg Fg
  for { set i 1 } { $i <= $Ns } { incr i } {
    set s state$i
    eval .c create oval \
      [expr $Cx($s)-$Vx($s)] [expr $Cy($s)-$Vy($s)] \
      [expr $Cx($s)+$Vx($s)] [expr $Cy($s)+$Vy($s)] \
      -width 0 -fill $Bg($s) -outline $Bg($s) -tags \{ $s state \}
    set tag to$s
    .c create line $Sx($s) $Sy $Cx($s) $Cy($s) -tags $tag \
      -fill $Fg($s)
  }
}
drawInitStates
#------------------------------------------------------------------------------
# stochastic random data generation from the model
#------------------------------------------------------------------------------
proc genPoints { } {
  global Np PointList Ns Trans Vx Vy Cx Cy Fg Trans
  eval .c delete [.c find withtag data]
  eval .c delete [.c find withtag line]
  set PointList {}
  set Np 0
  set prob 0
  for { set k 1 } { $k <= $Ns } { } {
    set i state$k
#    set t [expr sqrt(-2.*log(1.-[rand]))]

    set t [expr \
  [rand]+[rand]+[rand]+[rand]+[rand]+[rand]+[rand]+[rand]+[rand]+[rand]+[rand]+[rand]-6.]

    set u [expr 2*3.14159*[rand]]
    set x [expr $Cx($i)+$Vx($i)*$t*cos($u)]
    set y [expr $Cy($i)+$Vy($i)*$t*sin($u)]
    addPoint $x $y $Fg($i)
    update idletasks
    set prob [expr $prob+[logGauss $x $y $Cx($i) $Cy($i) $Vx($i) $Vy($i)]]
    if { [rand] < $Trans($i) } {
      set prob [expr $prob+log($Trans($i))]
      incr k
    } else {
      set prob [expr $prob+log(1.-$Trans($i))]
    }
  }
  drawShape
#  puts $prob
  .f.score config -text "log prob = $prob"
  barScore $prob
  return $prob
}
#------------------------------------------------------------------------------
# Viterbi path search
#------------------------------------------------------------------------------
proc Viterbi { } {
  global Ns Np PointList Cx Cy Vx Vy Trans Fg
  ### initialize the trellis ###
  for { set i 1 } { $i <= $Ns } { incr i } {
    set j state$i
    set score($j) -100000000.
  }
  set ts state1/0
  set back($ts) 0
  ### recursive computation ###
  set time 0
  foreach id $PointList {
    incr time
    set tags [lindex [.c itemconfigure $id -tags] 4]
    set x [lindex $tags 0]
    set y [lindex $tags 1]
#    puts stdout "$id (x, y) = $x, $y"
    for { set i $Ns } { $i > 1 } { incr i -1 } {
      set curr state$i
      set prev state[expr $i-1]
      set t [expr $score($prev)+log($Trans($prev))]
      set s [expr $score($curr)+log(1.-$Trans($curr))]
      set p [logGauss $x $y $Cx($curr) $Cy($curr) $Vx($curr) $Vy($curr)]
      if { $t > $s } {
        set score($curr) [expr $t+$p]
        set ts $curr/$time
        set back($ts) $prev
      } else {
        set score($curr) [expr $s+$p]
        set ts $curr/$time
        set back($ts) $curr
      }
    }
    if { $time == 1 } {
      set score(state1) [expr [logGauss $x $y \
        $Cx(state1) $Cy(state1) $Vx(state1) $Vy(state1)]]
    } else {
      set score(state1) [expr $score(state1)+[logGauss $x $y \
       $Cx(state1) $Cy(state1) $Vx(state1) $Vy(state1)]+log(1.-$Trans(state1))]
    }
    set ts state1/$time
    set back($ts) state1
  }
  set k state$Ns
  set score($k) [expr $score($k)+log($Trans(state1))]
#  puts $score($k)
  .f.score config -text "log prob = $score($k)"
  barScore $score($k)

#  ### show back pointers ###
#  for { set i 1 } { $i <= $time } { incr i } {
#    for { set j 1 } { $j <= $Ns } { incr j } {
#      set k state$j/$i
#      puts stdout back($k)=$back($k)
#    }
#  }

  ### Viterbi trace back ###
  set s state$Ns
  for { set t $time } { $t >= 1 } { incr t -1 } {
    .c itemconfigure [lindex $PointList [expr $t-1]] -fill $Fg($s)
    set st $s/$t
    set s $back($st)
#    puts stdout "back to $s/($t-1)"
  }
}

### display the Viterbi score ###
proc initbarScore { } {
  global W H
  .c create line $W $H $W $H -width 10 -tag score -fill #eb5808
}
proc barScore { score } {
  global W H Ns
  .c coords [.c find withtag score] $W $H $W [expr -$score*$H/$Ns/80]
}
initbarScore
#------------------------------------------------------------------------------
# logarithmic Gaussian probability density
#------------------------------------------------------------------------------
proc logGauss { x y xmean ymean xvar yvar } {
  return [expr \
  -pow(($x-$xmean)/$xvar,2.)/2.-pow(($y-$ymean)/$yvar,2.)/2.-log(abs($xvar*$yvar)) \
  ]
}
#------------------------------------------------------------------------------
# random data generation
#------------------------------------------------------------------------------
### Wichmann and Hill's random number generation algorithm ###
set Randx 1.
set Randy 1.
set Randz 1.
proc rand {} {
  global Randx Randy Randz
  set Randx [expr 171.*fmod($Randx,177.)-2.*floor($Randx/177.)]
  set Randy [expr 172.*fmod($Randy,176.)-35.*floor($Randy/176.)]
  set Randz [expr 170.*fmod($Randz,178.)-63.*floor($Randz/178.)]
  if { $Randx < 0 } { set Randx [expr $Randx+30269.] }
  if { $Randy < 0 } { set Randy [expr $Randy+30307.] }
  if { $Randz < 0 } { set Randz [expr $Randz+30323.] }
  set r [expr $Randx/30269.0+$Randy/30307.0+$Randz/30323.0]
  while { $r >= 1. } { set r [expr $r-1] }
#  if { $r <= 0. || $r >= 1. } { puts "Rand Error : rand = $r" }
  return $r
}
### randomize the seed with the current time (min and sec) ###
set Randx [exec /bin/date +%M%S]
#------------------------------------------------------------------------------
# end of the "wish" (tcl/tk) script: ShapeDesign
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
#  programming notes
#------------------------------------------------------------------------------
#
#   Points ... list of squares, each square containing a tag list
#              { x y "data" "move" }
#   
#
#
#