
package DAFS_Run;

use Exporter;
@ISA = (Exporter);
@EXPORT = qw(dafs_run);

sub new{

    my $this = shift;

    my (
	$tmpdir,
	$cosgap_model
	) = @_;

    my $DAFS_Run = {
	"tmp" => $tmpdir,
	"cosgap_model" => $cosgap_model,
    };

    bless $DAFS_Run, $this;
    return $DAFS_Run;
}

sub dafs_run
{
    my ($this, $rList, $rSimMatrix, $rgeneNamesIndex, $options) = @_;
    my $count=0;
    for(my $i=0; $i<=$#{$rList}; $i++)
    {
	$$rgeneNamesIndex[$i]=$$rList[$i][0];
	for(my $j=0; $j<=$i; $j++)
	{
	    if($i == $j){($$rSimMatrix[$i][$j][0], $$rSimMatrix[$i][$j][0]) = (0,0); next;}
	    my $tempfafile="$this->{tmp}/seq.fa";
	    my $tempcovfile="$this->{tmp}/seq.cov";
	    $this->makeSeqAndCovFile($i,$j,$rList,$tempfafile,$tempcovfile);
	    my $resultfile="$this->{tmp}/result.txt";
	    system("dafs $options $tempfafile $tempcovfile > $resultfile");
	    system("cat $resultfile >> $this->{tmp}/result2.txt");

	    @{$$rSimMatrix[$i][$j]} = $this->returnSimilarities($resultfile);
	    @{$$rSimMatrix[$j][$i]} = @{$$rSimMatrix[$i][$j]};
	    $count++;
	    if($count % 1000 == 0){print "$count pairs done\n";}
	}
    }
}

sub makeSeqAndCovFile
{
    my ($this, $i,$j,$rList,$tempfafile,$tempcovfile) = @_;
    open(OUT,"> $tempfafile") || die "Can't open : $tempfafile $!";
    print OUT ">$$rList[$i][0]\n"; 
    print OUT "$$rList[$i][1]\n"; 
    print OUT ">$$rList[$j][0]\n";
    print OUT "$$rList[$j][1]\n";
    close(OUT);

    #print "test : $$rList[$i][0] and $$rList[$j][0]\n"; 
  
    open(OUT,"> $tempcovfile") || die "Can't open : $tempcovfile $!";
    print OUT "$$rList[$i][2]\n"; 
    print OUT "$$rList[$j][2]\n";
    close(OUT);
}

sub returnSimilarities
{
    my ($this, $resultfile) = @_;
    open(RESULT,"$resultfile") || die "Can't open : $resultfile $!";

    my $temp = <RESULT>; #ignore first line of DAFS output
    my @SimAlign = $this->returnSimilarityOfAlign();
    $temp = <RESULT>; #ignore secondary structure alignment line of DAFS output
    my $SimShape = $this->returnMSsimilarity();

    push(@SimAlign, $SimShape);
    return (@SimAlign);
}

sub returnSimilarityOfAlign
{ 
    my $this = shift;
    my $temp = <RESULT>; chomp($temp); # ex. $temp = ">SS_cons 0.0110586 0.274345 .." 
    my @temp2=split(/\s+/,$temp);
    return (@temp2[1..$#temp2]);
}

sub returnMSsimilarity
{
    my $this = shift;

    my $firstGeneName=<RESULT>;
    my $firstGeneSeqAlignment=<RESULT>;
    my $firstGeneCov=<RESULT>; chomp($firstGeneCov);

    my $secondGeneName=<RESULT>;
    my $secondGeneSeqAlignment=<RESULT>;
    my $secondGeneCov=<RESULT>; chomp($secondGeneCov);

    if($this->{cosgap_model} eq "gap_free")
    {
	return $this->makeMSsimilarity2($firstGeneCov, $secondGeneCov);
    }
    elsif($this->{cosgap_model} eq "gap_fill")
    {
	return $this->makeMSsimilarity($firstGeneCov, $secondGeneCov);
    }
}

sub makeMSsimilarity
{
    my ($this, $firstGeneCov, $secondGeneCov)=@_;
    my @firstGeneCovs=split(/,/,$firstGeneCov);
    my @secondGeneCovs=split(/,/,$secondGeneCov);

    $this->CalGapCoverages(\@firstGeneCovs);
    $this->CalGapCoverages(\@secondGeneCovs);

    return $this->returnCosineCoeffient(\@firstGeneCovs,\@secondGeneCovs);
}

sub CalGapCoverages
{
    my ($this, $GeneCovs) = @_;
    my @GapStartPos=();
    my @GapEndPos=();

    $this->makeGapStartEndPos($GeneCovs, \@GapStartPos, \@GapEndPos);
    $this->putGapValueInCovs($GeneCovs, \@GapStartPos, \@GapEndPos); 
}

sub makeGapStartEndPos
{
    my ($this, $GeneCovs, $rGapStartPos, $rGapEndPos) = @_;
    my $FlagGapStart=0;
    for(my $i=0; $i<=$#{$GeneCovs}; $i++)
    {
	if($$GeneCovs[$i] eq "-" && $FlagGapStart == 0){ push(@$rGapStartPos,$i); $FlagGapStart = 1; }
	elsif($$GeneCovs[$i] ne "-" && $FlagGapStart == 1){ push(@$rGapEndPos,$i-1); $FlagGapStart = 0; }
	
	if($i == $#{$GeneCovs} && $FlagGapStart == 1){ push(@$rGapEndPos,$i); }
    }
}

sub putGapValueInCovs
{
    my ($this, $GeneCovs, $rGapStartPos, $rGapEndPos) = @_; 
    for(my $i=0; $i<=$#{$rGapStartPos}; $i++)
    {
	if($$rGapStartPos[$i] == 0)
	{
	    for(my $j=0; $j<=$$rGapEndPos[$i]; $j++)
	    {
		if($bothendzero){ $$GeneCovs[$j]=0;}
		else{ $$GeneCovs[$j]=$$GeneCovs[$$rGapEndPos[$i]+1]; }
	    }
	}
	elsif($$rGapEndPos[$i] == $#{$GeneCovs})
	{
	    for(my $j=$$rGapStartPos[$i]; $j<=$#{$GeneCovs}; $j++)
	    {
		if($bothendzero){ $$GeneCovs[$j]=0;}
		else{ $$GeneCovs[$j]=$$GeneCovs[$$rGapStartPos[$i]-1]; }
	    }
	}
	else
	{
	    for(my $j=$$rGapStartPos[$i]; $j<=$$rGapEndPos[$i]; $j++)
	    {
		$$GeneCovs[$j]=( $$GeneCovs[$$rGapStartPos[$i]-1]+$$GeneCovs[$$rGapEndPos[$i]+1] ) / 2;
	    }
	}
    }
}

sub makeMSsimilarity2
{
    my ($this, $firstGeneCov, $secondGeneCov)=@_;
    my @firstGeneCovs=split(/,/,$firstGeneCov);
    my @secondGeneCovs=split(/,/,$secondGeneCov);

    my @AlignPos = $this->getAlignPos(\@firstGeneCovs,\@secondGeneCovs);
    if($#AlignPos == -1){return 0;}

    @firstGeneCovs=@firstGeneCovs[@AlignPos];
    @secondGeneCovs=@secondGeneCovs[@AlignPos];
    return $this->returnCosineCoeffient(\@firstGeneCovs,\@secondGeneCovs);
}

sub getAlignPos
{
    my ($this, $rfirstGeneCovs, $rsecondGeneCovs) = @_;
    my @AlignPos=();
    for(my $i=0; $i<=$#{$rfirstGeneCovs}; $i++)
    {
	if($$rfirstGeneCovs[$i] ne "-" && $$rsecondGeneCovs[$i] ne "-"){push(@AlignPos,$i);}
    }
    return (@AlignPos);
}

sub returnCosineCoeffient
{
   my ($this, $rfirstGeneCovs, $rsecondGeneCovs) = @_;
   my $CosineCoeffient=0;

   my $sqrtsum_rfirstGeneCovs=$this->returnSqrtSum($rfirstGeneCovs);
   my $sqrtsum_rsecondGeneCovs=$this->returnSqrtSum($rsecondGeneCovs);
   my $Inner = $this->returnInner($rfirstGeneCovs,$rsecondGeneCovs);

   if($sqrtsum_rfirstGeneCovs == 0 || $sqrtsum_rsecondGeneCovs == 0){return 0;}
   return $CosineCoeffient= $Inner / ($sqrtsum_rfirstGeneCovs * $sqrtsum_rsecondGeneCovs);
}

sub returnSqrtSum
{
    my ($this, $rGeneCovs) = @_;
    my $sqrtsum=0;
    for(my $i=0; $i<=$#{$rGeneCovs}; $i++)
    {
	$sqrtsum += $$rGeneCovs[$i] ** 2;
    }
    return sqrt $sqrtsum;
}

sub returnInner
{
    my ($this, $rfirstGeneCovs,$rsecondGeneCovs) = @_;
    my $Inner=0;
    for(my $i=0; $i<=$#{$rfirstGeneCovs}; $i++)
    {
	$Inner+=$$rfirstGeneCovs[$i]*$$rsecondGeneCovs[$i];
    }
    return $Inner;
}

1;
