#!/usr/bin/perl

#Copyright (C) 2006-2008 Keio University
#(Kris Popendorf) <comp@bio.keio.ac.jp> (2006)
#
#This file is part of Murasaki.
#
#Murasaki 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 3 of the License, or
#(at your option) any later version.
#
#Murasaki 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 Murasaki.  If not, see <http://www.gnu.org/licenses/>.

###############
## simpal proof of concept perl implementation
###############

use strict;
use Getopt::Long;
use File::Basename;
use Pod::Usage;
use Data::Dumper;

BEGIN {
  unshift(@INC,(fileparse($0))[1].'perlmodules');
}
use Murasaki;

our $seedLength=3;
our $maxDistance=300;
our $tolerance=1;
our $maxLoop=300;

my ($help,$man);
GetOptions('help|?' => \$help, man => \$man, 'tolerance=i' => \$tolerance,
	   'maxLoop=i' =>\$maxLoop, 'seedLength=i' => \$seedLength,
	  'maxDistance=i' => \$maxLoop);
pod2usage(1) if $help or $#ARGV<1;
pod2usage(-exitstatus => 0, -verbose => 2) if $man;

our @seqs;
our $srcinit;

#find initial palidromes
foreach my $seq (@ARGV){
  my $start=time;
  die "File not found $seq" unless -e $seq;
  my ($fwd,$rev)=(lc `$root/geneparse -c -q $seq`,
		  lc `$root/antisense.pl -c $seq`);
  my %hash;
  fillHash($fwd,\%hash);
  fillHash($rev,\%hash,-1);
  my @pals=findPals(%hash);
  my %distmap=palsToDistmap(@pals);
  push(@seqs,{ pals => \@pals,
	       distmap => \%distmap,
	       file => $seq
       });
  $srcinit=0;
  print "Total palindromes: ".@pals."\n";
  print "Time to hash $seq: ".(time-$start)."\n";
}

my @srcvec=countvec($seqs[0]->{distmap});
my $src=$seqs[0];

foreach my $seq (@seqs[1..$#seqs]){
  my (@veca,@vecb);
  my ($namea,$nameb)=($src->{file},$seq->{file});
  my $similarity;
  my $start=time;
  print "Comparing to $namea to $nameb\n";
  if(0){
    print "Filtering matches to meet tolerance criteria (misses <= ".($tolerance<0 ? "any":$tolerance).")..\n";
    filterMatches($src,$seq,\@veca,\@vecb);
    print "Time spent filtering: ".(time-$start)."\n";
    $start=time;
    print "Vectorizing $namea...\n";
    @veca=countvec(@veca);
    print "Vectorizing $nameb...\n";
    @vecb=countvec(@vecb);
    print "Normalizing...\n";
    #  @veca=normalize(@veca);
    #  @vecb=normalize(@vecb);
    print "Dotting ...\n";
    $similarity=dot(\@veca,\@vecb);
  }else{
    $similarity=onepassCompute($src,$seq);
    print "$nameb matched $seq->{matched} times.\n";
  }
  print "Time spent: ".(time-$start)."\n";
  print "Similarity of $namea to $nameb: $similarity\n";
}

sub shortf(){
  my $a=shift;
  return sprintf("%.2f",$a);
}

sub onepassCompute(){
  my $total=0;
  my ($src,$target)=@_;
  foreach my $apal (@{$src->{pals}}){
    foreach my $bpal (@{$target->{pals}}){
      my $dist=abs($apal->{dist}-$bpal->{dist});
      next if $dist>$maxDistance;
#      print "Compare: ".join(" to ",prettyPal($apal),prettyPal($bpal))."\n";
      if($tolerance<0 or palMismatches($apal,$bpal)<=$tolerance){
#      if($tolerance<0 or mismatches($apal->{key},$bpal->{key})<=$tolerance){
	my $term=1.0/(exp($dist/5));
#	print "Term: $term\n";
	$total+=$term;
	$target->{matched}++;
      }
    }
  }
  return $total;
}

sub filterMatches(){
  my (%a,%b);
  my ($src,$target,$aout,$bout)=@_;
  foreach my $dist (keys(%{$src->{distmap}})){
    foreach my $apal (@{$src->{distmap}->{$dist}}){
      foreach my $bpal (@{$target->{distmap}->{$dist}}){
	if($tolerance<0 or mismatches($apal->{key},$bpal->{key})<=$tolerance){
	  $a{$apal}=$apal;
	  $b{$bpal}=$bpal;
	}
      }
    }
  }
  push(@$aout,values(%a));
  push(@$bout,values(%b));
}

sub dot { #compute dot product
  my ($a,$b)=@_;
  my $sum=0;
  for(0..($#$a<$#$b ? $#$a:$#$b)){
    $sum+=$$a[$_]*$$b[$_];
  }
  return $sum;
}

sub magnitude {
  my $sum;
  foreach(@_){
    $sum+=$_*$_;
  }
  return sqrt($sum);
}

sub vecDiv {
  my $div=shift;
  return map {$_/$div} @_;
}

sub normalize {
  return vecDiv(magnitude(@_),@_);
}

sub countvec {
  my @ret;
  foreach(@_){
    $ret[$_->{dist}]++;
  }
  return @ret;
}

sub max {
  my $best=pop;
  foreach(@_){
    $best=$_ if $best<$_;
  }
  return $best;
}

sub min {
  my $best=pop;
  foreach(@_){
    $best=$_ if $best>$_;
  }
  return $best;
}

sub palsToDistmap {
  my %distmap;
  foreach(@_){
    push(@{$distmap{$_->{dist}}},$_);
  }
  return %distmap;
}

sub palMismatches{
  my ($a,$b)=@_;
  return min(mismatches($a->{key},$b->{key}),
	     mismatches(revcomp($a->{key}),revcomp($b->{key})));
}

sub mismatches{
  my @a=split(//,pop);
  my @b=split(//,pop);
  my $count=@a;
  $count=@b if @b>@a; 
  for my $i (0..($a<$b ? $#a:$#b)){
    $count-- if $a[$i] eq $b[$i];
  }
  return $count;
}

sub findPals {
  my %in=@_;
  my @out;
  foreach my $key (keys %in){
    my @hits=@{$in{$key}};
    foreach my $rev (grep {$_<0} @hits){
      foreach my $fwd (grep {$_>0} @hits){
	my $pal={left=>$rev,
		 right => $fwd,
		 dist => $fwd+$rev,
		 key => $key};
	push(@out,$pal) if $pal->{dist}<=$maxLoop and $maxLoop>=0;
      }
    }
  }
  return @out;
}

sub prettyPal {
  my $pal=shift;
  return $pal->{left}." ".revcomp($pal->{key})." -{".$pal->{dist}."}- ".$pal->{key}." ".$pal->{right};
}

sub revcomp {
  ($_)=@_;
  y/agtcAGTC/tcagTCAG/;
  return reverse($_);
}

sub fillHash {
  my $seq=shift;
  my $hashref=shift;
  my %hash=%{$hashref};
  my $mult=(shift() ? -1:1);
  for(0..(length($seq)-$seedLength-1)){
    push(@{$hash{substr($seq,$_,$seedLength)}},($_+1)*$mult);
  }
  %{$hashref}=%hash;
}

__END__

=head1 NAME

simpal -- proof of concept

=head1 SYNOPSIS

simpal.pl [options] <ref seq> <target seq> [target2 target3 ...]

=head1 OPTIONS

 Set options by something like --tolerance=5 or -t 5 or whatever.
maxLoop -- specifies max separation for the ends of the palindromes
maxDistance -- specifies the maximum difference in palindrome length
               when comparing palindromes
tolerance -- max permissable misses in palindromes between sequences

Setting either to -1 means no limit

seedlength -- length of seeds to use for finding palindromes

=cut

