// sequence alignment

// X[1..M]  Y[1..N]  OPT[0..M, 0..N]
// del = gap penalty,  alpha = substitution penalty

// example X, Y are strings of English alphabets (lower case)
// alpha is 0 when X[i] Y[j] matched
//          1 when both are the same type (vowel or consonant)
//          3 when both are different type
// del = 2

del = 2
M = 4
N = 4

to prints s | i =
  i = 1
  while s[i] != 0 
    printc s[i]
    i = i + 1

// convert 2D to 1D
to index i j = (i*(M+1))+j

to prOPT | i j =
  for i 0 M 
    for j 0 N
      print OPT[index i j] space
    nl

// vowel is a, e, i, o, u
to isVowel c = 
  (c == 97) | (c == 101) | (c == 105) | (c == 111) | (c == 117)

// difference of the same type cost 1
// difference of the different type cost 3
to alpha i j | x y = 
  x = isVowel X[i]
  y = isVowel Y[j]  
  if X[i] == Y[j] 
    0 
  else if and (x == 1) (y == 1) 
    1
  else if and (x == 0) (y == 0)
    1
  else
    3

// "mean" = { 109, 101, 97, 110 }
// "name" = { 110, 97, 109, 101 }
to init = 
  X = array (M+2)
  Y = array (N+2)
  OPT = array ((M+1)*(N+1))
  X[0] = 0
  X[1] = 109
  X[2] = 101
  X[3] = 97
  X[4] = 110
  X[5] = 0
  Y[0] = 0
  Y[1] = 110
  Y[2] = 97
  Y[3] = 109
  Y[4] = 101
  Y[5] = 0

to min2 a b = if a < b a else b
to min a b c = min2 (min2 a b) (min2 b c)

to align | i j =
  for i 1 M
    OPT[index 0 i] = i*del
  for j 1 N
    OPT[index j 0] = j*del
  for i 1 M
    for j 1 N
      OPT[index i j] = 
      min ((alpha i j) + OPT[index (i-1) (j-1)])
      (del + OPT[index (i-1) j])
      (del + OPT[index i (j-1)])
 
to trace i j | x y d v =
  v = OPT[index i j]
  print i space
  print j nl
  if or (i==0) (j==0) break
  d = (alpha i j) + OPT[index (i-1) (j-1)]
  x = del + OPT[index (i-1) j]
  y = del + OPT[index i (j-1)]
  if v == d       trace (i-1) (j-1)
  else if v == x  trace (i-1) j
  else if v == y  trace i (j-1)

// space efficient alignment 
// B[0..M, 0..1] collapse column of OPT[.] 
// B[i,0] = OPT[i,j-1]   B[i,1] = OPT[i,j]

to align2 | i j =
  B = array ((M+1)*2)
  for i 1 M
    B[index i 0] = i*del
  for j 1 N
    B[index 0 j] = j*del
    for i 1 M
      B[index i j] = 
      min ((alpha i j) + B[index (i-1) 0])
      (del + B[index (i-1) 1])
      (del + B[index i 0])
    for i 1 M
      B[index i 0] = B[index i 1]

to main = 
  init
  prints X space prints Y nl
//  align
//  prOPT
//  trace M N
  align2
  print B[index M 0] nl

main