#!/usr/bin/perl
use strict;
use warnings;

#debug = 1 displays score and traceback matrix
my $debug = 0;
my $out_fh=\*STDOUT;
my @seq1;
my @seq2;

# Nutzung:
# Scriptname [debug] <Seq_File1> <Seq_File2> <Ausgabe_Datei>
#
# wenn "debug" dann Anzeige von score und traceback Matrix
#
# wenn "<Seq_File1>" oder "<Seq_File2>" gleich "-" werden sie ignoriert
# und nach der Sequenz gefragt
#
# wenn "<Ausgabe_Datei>" nicht Angegeben wird auf STDOUT geschrieben

if(@ARGV)
{
  if($ARGV[0] eq 'debug')
  {
    $debug = 1;
    shift(@ARGV);
  }

  @seq1=read_seq(shift(@ARGV));
  @seq2=read_seq(shift(@ARGV));

  if(@ARGV)
  {
    my $out_file=shift(@ARGV);
    if($out_file && open(my $fh, '>', $out_file))
    { $out_fh=$fh; }
    else
    { warn("Kann Augabedatei $out_file nicht öffnen! ($!)\n"); }
  }
}


@seq1=ask_seq("Gib erste Sequenenz ein") unless(@seq1);
@seq2=ask_seq("Gib zweite Sequenenz ein") unless(@seq2);

my $now = time;

####################################################################################################
#                      Schema für Score                                                            #
####################################################################################################

my $MATCH    =  1; # +1 für Buchtstabe, der ein Match macht
my $MISMATCH = -1; # -1 für Buchstabe, der einen Mismatch macht
my $GAP      = -1; # -1 für Lücke


####################################################################################################
#                      Initialisierung / Programmstart                                             #
####################################################################################################
my @matrix;
$matrix[0][0]{score}   = 0;
$matrix[0][0]{pointer} = "none";
for(1..@seq1)
{
  $matrix[0][$_]{score}   = $GAP * $_;
  $matrix[0][$_]{pointer} = "left";
}
for(1..@seq2)
{
  $matrix[$_][0]{score}   = $GAP * $_;
  $matrix[$_][0]{pointer} = "up";
}

####################################################################################################
#                                   Auffüllen                                                      #
####################################################################################################
for my $i (1..@seq2)
{
  for my $j (1..@seq1)
  {
    ####################################################################################################
    #                            berechnet match score                                                 #
    ####################################################################################################
    my $diagonal_score = $matrix[$i-1][$j-1]{score};
    $diagonal_score += $seq1[$j-1] eq $seq2[$i-1]? $MATCH: $MISMATCH;


    ####################################################################################################
    #                           berechnet gap scores                                                   #
    ####################################################################################################
    my $up_score   = $matrix[$i-1][$j]{score} + $GAP;
    my $left_score = $matrix[$i][$j-1]{score} + $GAP;



    ####################################################################################################
    #                             aussuchen des best score                                             #
    ####################################################################################################
    if ($diagonal_score >= $up_score)
    {
      if ($diagonal_score >= $left_score)
      {
        $matrix[$i][$j]{score}   = $diagonal_score;
        $matrix[$i][$j]{pointer} = "diagonal";
      }
      else
      {
        $matrix[$i][$j]{score}   = $left_score;
        $matrix[$i][$j]{pointer} = "left";
      }
    }
    else
    {
      if ($up_score >= $left_score)
      {
        $matrix[$i][$j]{score}   = $up_score;
        $matrix[$i][$j]{pointer} = "up";
      }
      else
      {
        $matrix[$i][$j]{score}   = $left_score;
        $matrix[$i][$j]{pointer} = "left";
      }
    }
  }
}

##########################################################################################################
#              print matrix                                                                              #
##########################################################################################################

if ($debug ==1 )
{
  print "The Score Matrix\n";
  print "-----------------------\n";
  print "\t";
  print $_."\t" for (@seq1);
  print "\n";
  for my $u (0..@seq2)
  {
    print $u==0?"*":$seq2[$u-1] ,"\t";
    print $matrix[$u][$_]{score}."\t" for(0..@seq1);
    print "\n";
  }

  #####################################################################################################
  print "The Trace Matrix\n";
  print "-----------------------\n";
  print "\t";
  print $_."\t" for (@seq1);
  print "\n";

  for my $u (0..@seq2)
  {
    print $u==0?"*":$seq2[$u-1] ,"\t";
    for my $v (0..@seq1)
    {
      my $x=$matrix[$u][$v]{pointer};
      print  "N\t" if ($x eq "none");
      print  "L\t" if ($x eq "left");
      print  "U\t" if ($x eq "up");
      print  "D\t" if ($x eq "diagonal");
    }
    print "\n";
  }
}

my $align1 = "";
my $align2 = "";

my $i=@seq2;
my $j=@seq1;
my $score = $matrix[$i][$j]{score};
my $x=$matrix[$i][$j]{pointer};
while ($x ne "none" && $i>=0 && $j>=0)
{
  if ($x eq "diagonal")
  {
    $align1 .= $seq1[$j-1];
    $align2 .= $seq2[$i-1];
    $i--;
    $j--;
  }
  elsif ($x eq "left")
  {
    $align1 .= $seq1[$j-1];
    $align2 .= "-";
    $j--;
  }
  elsif ($x eq "up")
  {
    $align1 .= "-";
    $align2 .= $seq2[$i-1];
    $i--;
  }
  $x=$matrix[$i][$j]{pointer};
}

$align1 = reverse $align1;
$align2 = reverse $align2;

print $out_fh "Best Alignment\n";
print $out_fh "--------------\n";
print $out_fh "Seq1:\t$align1\n";
print $out_fh "Seq2:\t$align2\n";
print $out_fh "score = $score\n";

print time - $now." second\(s\)\n";

########################################################################
########################################################################
sub read_seq
{
  my $file=shift;
  return () unless($file);
  return () if($file eq '-');
  if(open( my $fh, '<', $file))
  {
    local $/=undef;
    chomp(my $str=<$fh>);
    close($fh);
    return split(//, $str);
  }

  warn("Kann Datei $file nicht öffnen! ($!)\n");
  return ();
}

sub ask_seq
{
  my $txt=shift // '';
  print "$txt\n";
  chomp(my $seq=<STDIN>);
  return split(//,$seq);
}