#+     Copyright (C) A.D. 2002 by Kevin C. O'Kane  
#+
#+     Kevin C. O'Kane
#+     Department of Computer Science
#+     University of Northern Iowa
#+     Cedar Falls, IA 50613-0507 USA
#+
#+     okane@cs.uni.edu
#+     anamfianna@earthlink.net
#+     http://www.cs.uni.edu/~okane
#+
#+
#+ This program is free software; you can redistribute it and/or modify
#+ it under the terms of the GNU General Public License as published by
#+ the Free Software Foundation; either version 2 of the License, or
#+ (at your option) any later version.
#+ 
#+ This program is distributed in the hope that it will be useful,
#+ but WITHOUT ANY WARRANTY; without even the implied warranty of
#+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#+ GNU General Public License for more details.
#+ 
#+ You should have received a copy of the GNU General Public License
#+ along with this program; if not, write to the Free Software
#+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
#+
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

^^Score(a,b)
      if a="-" quit "-1"
      if b="-" quit "-1"
      if a=b quit "2"
      quit "-1"

^^Cost(s,t,i,j)
      if $d(^a(i,j)) quit ^a(i,j)
      if i=0,j=0 set ^a(i,j)=0 quit "0"
      if i'=0,j=0 set x=$$^^Cost(s,t,i-1,0)+$$^^Score($e(s,i,i),"-") set ^a(i,j)=x q x
      if i=0,j'=0 set x=$$^^Cost(s,t,0,j-1)+$$^^Score("-",$e(t,j,j)) set ^a(i,j)=x q x
      set a=$$^^Cost(s,t,i-1,j-1)+$$^^Score($e(s,i,i),$e(t,j,j))
      set b=$$^^Cost(s,t,i-1,j)+$$^^Score($e(s,i,i),"-")
      set c=$$^^Cost(s,t,i,j-1)+$$^^Score("-",$e(t,j,j))
      set x=a
      if b>x set x=b
      if c>x set x=c
      set ^a(i,j)=x
      quit x

      zmain
      set s="acbcdb"
      set t="cadbd"
      kill ^a
      write !,"Smith-Waterman Optimal Sequence Aligments between ",s," and ",t,!!
      do ^^Cost(s,t,$l(s),$l(t))
      for i=0:1:$l(t) write ?(i+1*6),$j(i,5)
      write !," "
      for i=1:1:$l(t) write ?(i+1*6),$j($e(t,i,i),5)
      write !
      set i=-1
      for  do
      . set i=$n(^a(i))
      . if i="-1" break
      . write $j(i,3)," "
      . write:i>0 $e(s,i,i),":"
      . set j=-1
      . for  do
      .. set j=$n(^a(i,j))
      .. if j="-1" break
      .. write ?(j+1*6),$j(^a(i,j),5)
      . write !
      halt