#!/usr/bin/perl

$command_dir = "/home6/jpei/speltest/SPEL/bin";

$cpunum = 5;

if($#ARGV<0) {
	print "Usage: SPEL.perl alignment_file [-cpu n] \n";
	print "alignment_file: plain alignment or clustal format\n";
	exit(1);
}

#if($ARGV[0] =~ /LINUX/i) { $platform = "LINUX_RH_PC"; }
#elsif($ARGV[0] =~ /SOLARIS/i) { $platform = "SOLARIS_INTEL"; }
#elsif($ARGV[0] =~ /SPARC/i) { $platform = "SOLARIS_SPARC"; }
#else { print "Can not recognize your platform\n"; exit(1); }

$alnFile = $ARGV[0];

for($i=0;$i<$#ARGV;$i++) {
	if($ARGV[$i] eq "-cpu") {$cpunum = $ARGV[$i+1];}
}

$command = "$command_dir/ancescon"." -RO -C -i ".$alnFile." -o junk";
my $status = system($command);
if($status != 0) {
	print "Something is wrong when we construct the phylogenetic tree\n";
	exit(1);
}

$treeFile = $alnFile.".tre";

$command = "mpirun -np $cpunum $command_dir/simulation"." -i ".$alnFile." -t ".$treeFile." -o ".$alnFile.".sim";
print $command, "\n";
my $status = system($command);
if($status != 0) {
        print "Something is wrong when we do simulation\n";
        exit(1);
}

$simFile = $alnFile.".sim";

# read simulation output file
$flag = 0;
@Z = ();
@sim_Z = ();
open(COG, $simFile) or die "can not open file: $simFile to read\n";
while(<COG>) {
        chomp;
        if($_ =~ /site   entropy     avgEntropy  logL  avgLogL/) {
                $flag = 1;
                next;
        }
        if($_ =~ /Log likelihood for each site after shuffle/) {
                $flag = 2;
                next;
        }
        if($_ =~ /check/) {
                $flag = 0;
                last;
        }
        if($flag == 1) {
                push @Z, [split];
        }
        elsif($flag ==2) {
                push @sim_Z, $_;
        }
}
close(COG);

# record the total sites
$total_sites = $#Z+1;

# output file *.Z
$Zname = $alnFile.".Z";
open(COG, ">$Zname") or die "can not open file: $Zname to write\n";
select(COG);
for($i=0; $i<=$#Z; $i++) {
        printf("%10.4f\n", $Z[$i][3]);
}
close(COG);
select(STDOUT);
# output file *_sim.Z
$simZname = $alnFile."_sim.Z";
open(COG, ">$simZname") or die "can not open file: $simZname to write\n";
select(COG);
for($i=0; $i<=$#sim_Z; $i++) {
        print "$sim_Z[$i]\n";
}
close(COG);
select(STDOUT);

# Calculate p-value for each position
$command = "$command_dir/specificity"." -R -O ".$Zname." -M ";
$command = $command.$simZname." > ".$alnFile.".result";
system($command);

# read result file and display the result
$flag = 0;
@result = ();
$resultFile = $alnFile.".result";
open(COG, $resultFile) or die "can not open file $resultFile to read\n";
while(<COG>) {
	chomp;
	if($_ =~ /site    prob1  Likelihood/) {
		$flag = 1;
		next;
	}
	if($flag == 1) {
		if($_ =~ /^$/) {last;}
		
		push @result, [split];
	}
}
close(COG);

$specificityFile = $alnFile.".specificity";
open(COG, "> $specificityFile") or die "can not open file $specificityFile to write\n";
select(COG);
printf("Rank  Site    P-Value\n");
for $i (0..$#result) {
	printf("%4d   %4d   %12.6e\n",$i, $result[$i][0],$result[$i][1]);
}
close(COG);
select(STDOUT);
print "Thank you for using SPEL program!\n";
