
use strict;
my $file1 = $ARGV[0];
my $file2 = $ARGV[1];
my ($ref1, $ref2, $W1, $W2);

#read first motif from a simple text file
($ref1, $W1) = readMotifFromSimpleFile($file1);
my @matrix1 = @$ref1;
#printMatrix($ref1, $W1);

#read second motif from a simple text file
($ref2, $W2) = readMotifFromSimpleFile($file2);
my @matrix2 = @$ref2;
#printMatrix($ref2, $W2);

my $dist = computeDistance($ref1, $W1, $ref2, $W2);
print sprintf("%.4f\n",$dist);

# compute the distance between two PSSMs
sub computeDistance()
{
   # read the matrices and make sure the first matrix is shoter
   my ($ref1, $W1, $ref2, $W2) = @_;
   my $reversed = 0;
   if ($W1 > $W2) {
      my $tmp = $W1; $W1 = $W2; $W2 = $tmp;
      $tmp = $ref1; $ref1 = $ref2; $ref2 = $tmp;
      $reversed = 1;
   }
   my @m1 = @$ref1;
   my @m2 = @$ref2;

   # set the constants
   my $MINOVERLAP = 3;
   my $REGOVERLAP = 6;
   my $ENTROPYCUTOFF = 1;
   
   # check the sizes and set the overlap as min($REGOVERLAP,$W1,$W2)
   if (scalar(@m1) != 4 || scalar(@m2) != 4) {
      print "computeDistance: One of the input matrices is not size 4!\n";
      return;
   }
   if ($W1 < $MINOVERLAP) {
      print "computeDistance: One of the input matrices is too small!\n";
      return; 
   }
   my $overlap = $REGOVERLAP;
   if ($W1 < $overlap) { $overlap = $W1; }
   
   # compute the reverse complement for the first matrix
   my @m1rev = reverse(@m1); # the lines are reversed
   for (my $i=0; $i<4; $i++) { # now reverse the columns
      my $tempref = $m1rev[$i];
      my @temparray = reverse(@$tempref);
      $m1rev[$i] = \@temparray;
   }
   
   # shift the matrices one against the other and
   # for every possible alignment compute the distance,
   # keeping track of the smallest distance
   my $dmin = 999999999999;
   my $len = $overlap;
   
   for ($len = $overlap; $len<=$W1; $len++) {
   for (my $i1=0; $i1<$W1-$len+1; $i1++) {
      for (my $i2=0; $i2<$W2-$len+1; $i2++) {
         # compute the distance between the matrices 
         # of size $len starting at $i1 and $i2
         my ($sum,$sumrev,$f,$frev,$E,$val) = (0,0,0,0,0,0);
         my $log2 = log(2);
         
         for (my $i=0; $i<$len; $i++) #for each column in the alignment
         {
            my ($s,$srev) = (0,0);
            for (my $k=0; $k<4; $k++) {#for each row (ACGT)
               $s += ($m1[$k][$i1+$i] - $m2[$k][$i2+$i])**2;
               $srev += ($m1rev[$k][$i1+$i] - $m2[$k][$i2+$i])**2;
               
               #we always compute the entropy for the first initial matrix
               if ($reversed == 1) { $val=$m2[$k][$i2+$i]; }
               else { $val=$m1[$k][$i1+$i]; } 
               if ($val==0) {
               	$val = 0.001;
               }
               $E += $val*log($val)/$log2;
            }
            if (sqrt($s/2) > 0.8) { $f++; } #the distance for any column
            if (sqrt($srev/2) > 0.8) { $frev++; } #should not be > 0.8
            $sum += sqrt($s);
            $sumrev += sqrt($srev);
         }
         $E = 2 + $E/$len;
         my $d = $sum / ($len*sqrt(2));
         #print "$f $len ", $f/$len, " ", 1/6, "\n";
         if (($f/$len >= 1/6) || ($E <= $ENTROPYCUTOFF)) { $d = 1; }
         if ($d < $dmin) {
            $dmin = $d;
         }
         
         $d = $sumrev / ($len*sqrt(2));
         #print "$f $len ", $f/$len, " ", 1/6, "\n";
         if (($frev/$len >= 1/6) || ($E <= $ENTROPYCUTOFF)) { $d = 1; }
         if ($d < $dmin) {
            $dmin = $d;
         }
      }#for $i2
   }#for $i1
  }
   
   return $dmin;
}

# read a motif from a simple text file (4 lines: acgt, W columns)
sub readMotifFromSimpleFile()
{
   my ($filename) = @_;
   my @matrix = ();
   my $W;
   open(IN, $filename) || die ("$filename does not exist\n");
   for (my $i=0; $i<4; $i++) {
     my $line = <IN>; chomp($line);
     my @array = split(" ", $line); $W = scalar(@array);
     push @matrix, \@array; 
   }
   close(IN);
   return (\@matrix, $W);
}

# print a matrix
sub printMatrix()
{
   my ($ref, $W) = @_;
   my @matrix = @$ref;
   foreach $ref (@matrix) {
      for (my $i=0; $i<$W; $i++) {
         print $ref->[$i], "* ";
      }
      print "*$W\n";
   }
}
