# IT Universitet i København - Algorithm Design Course (2010)
# ex5: Sequence Alignment
# Passes al the tests
# author: Alex Usbergo (aleu@itu.dk)

def fetch_data(matrix_filename, input_filename)
  f = File.open(matrix_filename)

  #reading the BLOSUM matrix
  blosum = Hash.new
  columns = f.readline.split[0..-1].map{|c| c.to_sym}

  while (!f.eof?) do
    values = f.readline.split[0..-1]
    row = values.shift.to_sym
    columns.each{|c| blosum[c] = {} if blosum[c].nil?; blosum[c][row] = values.shift.to_i}
  end

  f = File.open(input_filename)

  #reading the input file
  input = []
  while (!f.eof?) do
    input << {:name => f.readline.split[0][1..-1].to_sym,
              :data => ("*" + f.readline).split(//)[0..-2].map{|c| c.to_sym}}
  end

  return [blosum, input]
end

def seqalign(x, y, m, blosum)
  x.size.times{|i| m[i] = [-4*(i)]}
  y.size.times{|j| m[0][j] = -4*(j)}

  for i in 1..x.size-1 do
    for j in 1..y.size-1 do
      m[i][j] = [blosum[x[i]][y[j]] + m[i-1][j-1], -4 + m[i][j-1], -4 + m[i-1][j]].max
    end
  end
  m[x.size-1][y.size-1]
end

def seqprint(x, y, matrix, blosum)
  aligned_x = []
  aligned_y = []

  i, j = x.size-1, y.size-1

  while !i.zero? || !j.zero?

    #match
    if matrix[i][j] == matrix[i-1][j-1] + blosum[x[i]][y[j]]
      aligned_x << x[i] # =>
      aligned_y << y[j]
      j -= 1
      i -= 1

    #gaps
    elsif matrix[i][j] == -4 + matrix[i][j-1]
      aligned_x << :- # =>
      aligned_y << y[j]
      j -= 1

    elsif matrix[i][j] == -4 + matrix[i-1][j]
      aligned_x << x[i]  # =>
      aligned_y << :-
      i -= 1
    end
  end

  x_s = "" ; y_s = ""
  aligned_x.reverse_each{|c| x_s += c.to_s}
  aligned_y.reverse_each{|c| y_s += c.to_s}

  puts "#{x_s}\n#{y_s}\n\n"
end

#blosum, input = fetch_data("data/BLOSUM62.txt", "data/Toy_FASTAs.in")
blosum, input = fetch_data("data/BLOSUM62.txt","data/HbB_FASTAs.in")


for i in 0..input.size-1 do
  for j in i+1..input.size-1 do
    matrix = []
    value  = seqalign(input[i][:data], input[j][:data], matrix, blosum)
    puts "#{input[i][:name]}--#{input[j][:name]}: #{value}"
    seqprint(input[i][:data], input[j][:data], matrix, blosum)
  end
end