package BLAST::Util;
BEGIN{ push(@INC, $ENV{LOCUS_PERL_MODULES}); }
use Carp;
use warnings;
use strict;
use Utils::Array;
use FASTA::Util;
use FASTA::Sequence;
use Sequence::Locus;
use BLAST::Run;

###################################################
# Global variables
my $validExtensions = ['blast', 'bls'];

###################################################
sub getBlastFiles {
    my $dir = shift or croak "Directory name is not defined.";
    $dir .= '/' unless $dir =~ /\//;
    my $blastFiles = [];
    opendir DIR, $dir or die "Can't open $dir as a directory.";
    while( defined(my $file = readdir(DIR)) ){
	if( isValidBlastFile($file) ){
	    push(@$blastFiles, $dir.$file);
	}
    }
    closedir DIR;
    return $blastFiles;
}

###################################################
sub isValidBlastFile {
    my $filename = shift or croak "File name is not defined.";
    foreach my $extension (@$validExtensions){
	if( $filename =~ /\.$extension$/ ){
	    return 1; 
	    last;
	}
    }
    return 0;
}

###################################################
sub extractBestHits {
    my $hits = shift or croak "Blast hits are not defined.";
    my $queryIDToBestHit = {};
    foreach my $hit (@$hits){
	my $queryID = $hit->getQueryID();
	my $eValue = $hit->getEValue();
	if( $queryIDToBestHit->{$queryID} ){
	    if( $eValue < $queryIDToBestHit->{$queryID}->getEValue() ){
		$queryIDToBestHit->{$queryID} = $hit;
	    }
	} else {
	    $queryIDToBestHit->{$queryID} = $hit;
	}
    }
    my $bestHits = [];
    foreach my $queryID (keys %$queryIDToBestHit){
	my $bestHit = $queryIDToBestHit->{$queryID};
	push(@$bestHits, $bestHit);
    }
    return $bestHits;
}

###################################################
sub getQueryIDToHits {
    my $hits = shift or croak "Blast hits are not defined.";
    my $queryIDToHits = {};
    foreach my $hit (@$hits){
	my $queryID = $hit->getQueryID();
	if( $queryIDToHits->{$queryID} ){
	    push( @{$queryIDToHits->{$queryID}}, $hit );
	} else {
	    $queryIDToHits->{$queryID} = [$hit];
	}
    }
    return $queryIDToHits;
}

###################################################
sub extractReverseBestHits {
    my $hits = shift or croak "Blast hits are not defined.";
    my $subjectIDToBestHit = getSubjectIDToBestHit($hits);
    my $bestHits = [];
    foreach my $subjectID (keys %$subjectIDToBestHit){
	my $bestHit = $subjectIDToBestHit->{$subjectID};
	push(@$bestHits, $bestHit);
    }
    return $bestHits;
}

###################################################
sub getSubjectIDToHits {
    my $hits = shift or croak "Blast hits are not defined.";
    my $subjectIDToHits = {};
    foreach my $hit (@$hits){
	my $subjectID = $hit->getSubjectID();
	if( $subjectIDToHits->{$subjectID} ){
	    push( @{$subjectIDToHits->{$subjectID}}, $hit );
	} else {
	    $subjectIDToHits->{$subjectID} = [$hit];
	}
    }
    return $subjectIDToHits;
}

###################################################
sub getSubjectIDToBestHit {
    my $hits = shift or croak "Blast hits are not defined.";
    my $subjectIDToBestHit = {};
    foreach my $hit (@$hits){
	my $subjectID = $hit->getSubjectID();
	my $eValue = $hit->getEValue();
	if( $subjectIDToBestHit->{$subjectID} ){
	    if( $eValue < $subjectIDToBestHit->{$subjectID}->getEValue() ){
		$subjectIDToBestHit->{$subjectID} = $hit;
	    }
	} else {
	    $subjectIDToBestHit->{$subjectID} = $hit;
	}
    }
    return $subjectIDToBestHit;
}

###################################################
sub extractReverseNextBestHits {
    my $hits = shift or croak "Blast hits are not defined.";
    my $bestHits = shift or croak "Best hits are not defined.";
    my $nonBestHits = removeOverlappedHits($hits, $bestHits);
    return extractReverseBestHits($nonBestHits);
}

###################################################
sub removeOverlappedHits {
    my $hits1 = shift or croak "Blast hits (1) are not defined.";
    my $hits2 = shift or croak "Blast hits (2) are not defined.";
    my $remainedHits = [];
    foreach my $hit1 (@$hits1) {
	my $flgRemained = 1;
	foreach my $hit2 (@$hits2){
	    if( $hit1->isOverlappedOnSubject($hit2) ){
		$flgRemained = 0;
		last;
	    }
	}
	if( $flgRemained ){
	    push( @$remainedHits, $hit1 );
	}
    }
    return $remainedHits;
}

###################################################
sub extractLongHits {
    my $hits = shift or croak "Blast hits are not defined.";
    my $lengthCut = shift;
    my $longHits = [];
    foreach my $hit (@$hits){
	if( $hit->getAlignmentLength() >= $lengthCut ){
	    push(@$longHits, $hit);
	}
    }
    return $longHits;
}

###################################################
sub extractHighIdentityHits {
    my $hits = shift or croak "Blast hits are not defined.";
    my $cutoff = shift;
    my $hiHits = [];
    foreach my $hit (@$hits){
	if( $hit->getPercentIdentity() >= $cutoff ){
	    push( @$hiHits, $hit );
	}
    }
    return $hiHits;
}

###################################################
sub EvalueFilter {
    my $hits = shift or croak "Blast hits are not defined.";
    my $cutoff = shift;
    unless( $cutoff ){
	return $hits;
    }
    my $passedHits = [];
    foreach my $hit (@$hits){
	if( $hit->getEValue() <= $cutoff ){
	    push( @$passedHits, $hit );
	}
    }
    return $passedHits;
}

###################################################
sub getUniqueNrGIs {
    my $hits = shift or croak "Blast hits are not defined.";
    my $nrGIs = getNrGIs($hits);
    return Utils::Array::unique($nrGIs);
}

###################################################
sub getNrGIs {
    my $hits = shift or croak "Blast hits are not defined.";
    my $nrGIs = [];
    foreach my $hit (@$hits){
	my $nrGI = getNrGI($hit);
	push(@$nrGIs, $nrGI);
    }
    return $nrGIs;
}

###################################################
sub getNrGI {
    my $hit = shift;
    my $subjectID = $hit->getSubjectID();
    if( $subjectID =~ /gi.(\d+)./ ){
	return $1;
    }
    croak "No NrGI was found in the hit:\n"
	. $hit . "\n";
}

###################################################
sub getNonRedundantSeqs {
    my $hits = shift;
    my $seqs = shift;
    my $cutPercentIdentity = shift;
    my $cutPercentOverlap = shift;
    my $redundantTitles = getRedundantTitles
	( $hits, $seqs, $cutPercentIdentity, $cutPercentOverlap );
    my $nonRedundantSeqs = [];
    foreach my $seq (@$seqs){
	if( !defined($redundantTitles->{$seq->getTitle()}) ){
	    push(@$nonRedundantSeqs, $seq);
	}
    }
    return $nonRedundantSeqs;
}

###################################################
sub getRedundantTitles {
    my $hits = shift;
    my $seqs = shift;
    my $seqHash = FASTA::Util::indexing( $seqs );
    my $cutPercentIdentity = shift;
    my $cutPercentOverlap = shift;
    my $redundantTitles = {};
    foreach my $hit (@$hits){
	if(isRedundant
	   ($hit, $seqHash, $cutPercentIdentity, $cutPercentOverlap)){
	    my $redundantTitle = $hit->getQueryID();
	    if( length($seqHash->{$hit->getSubjectID()})
		< length($seqHash->{$hit->getQueryID()})){
		$redundantTitle = $hit->getSubjectID();
	    }
	    $redundantTitles->{$redundantTitle} = 1;
	}
    }
    return $redundantTitles;
}

###################################################
sub isRedundant {
    my $hit = shift;
    my $seqs = shift;
    my $cutPercentIdentity = shift;
    my $cutPercentOverlap = shift;
    if( isSame( $hit ) ){
	return 0;
    }
    if( isMatuallyIdentical( $hit, $cutPercentIdentity ) ){
	if( isMatuallyOverlap( $hit, $seqs, $cutPercentOverlap ) ){
	    return 1;
	}
    }
    return 0;
}

###################################################
sub isSame {
    my $hit = shift;
    if( $hit->getQueryID() eq $hit->getSubjectID() ){
	return 1;
    } else {
	return 0;
    }
}

###################################################
sub isMatuallyIdentical {
    my $hit = shift;
    my $cutPercentIdentity = shift;
    if( $hit->getPercentIdentity() > $cutPercentIdentity ){
	return 1;
    } else {
	return 0;
    }
}

###################################################
sub isMatuallyOverlap {
    my $hit = shift;
    my $seqs = shift;
    my $cutPercentOverlap = shift;
    my $percentOverlap = calcPercentOverlap( $hit, $seqs );
    if( $percentOverlap > $cutPercentOverlap ){
	return 1;
    } else {
	return 0;
    }
}

###################################################
sub calcPercentOverlap {
    my $hit = shift;
    my $seqs = shift;
    my $querySeqLength   = length($seqs->{$hit->getQueryID()});
    my $subjectSeqLength = length($seqs->{$hit->getSubjectID()});
    my $queryMatchLength = $hit->getQueryMatchLength();
    my $subjectMatchLength = $hit->getSubjectMatchLength();
    my $queryPercentOverlap = $queryMatchLength / $querySeqLength * 100;
    my $subjectPercentOverlap = $subjectMatchLength / $subjectSeqLength * 100;
    return Utils::Array::max($queryPercentOverlap, $subjectPercentOverlap);
}

###################################################
sub getReferenceLoci {
    my $hits = shift or croak "Blast hits are not defined.";
    my $loci = [];
    foreach my $hit (@$hits){
	my $locus = $hit->getSubjectLocus();
	push(@$loci, $locus);
    }
    return $loci;
}

###################################################
sub runBlast {
    my $params = shift;
    my $prefix = shift;
    my $job    = BLAST::Run->new();
    $job->setConfig( $params );
    $job->run( $prefix );
}

1;
