#!/usr/bin/perl -w

#################################################################################
# This script starts $NbMaxClusters jobs, varying the number of clusters used   #
# in the AHRS algorithm implemented in PAML (ver. 3.15)                         #
#                                                                               #
# Comments to: sarisbro at uottawa dot ca                                       #
# Date: May. 2006                                                               #
#################################################################################

use strict;
use File::Copy;
use Math::Random;
use LWP::Simple;
use Parallel::ForkManager;

my $NbMaxClusters = 2;
my $NbReps = 1000;
my $SGE = 1;

#################################################################################
# MAIN                                                                          #
#################################################################################

my($OStype,$pwd) = get_pwd();
my $curdir = ''; if($curdir eq ''){$curdir = $pwd;}else{$curdir = $ARGV[0];}
my $PI = 3.141592653;
my $MAX_PROCESSES = 2; # 4
my $pm = new Parallel::ForkManager($MAX_PROCESSES);
random_set_seed(time(),int(time()/$PI));
my($tree_file_name,$data_file_name,@tree_file_line,@seq_file_line,$ClusMethod,@ClusMethod_line);
my($ndata,@ndata_line,$subModel,@subModel_line,$start_at,$get_SE,$sge,$fix_kappa,$fix_alpha,$ncatG,$Small_Diff,$cleandata);
my($alpha,$kappa);

#if(!(-e "mybaseml")){
#	system "make mybaseml";
#}

if(!(-e "baseml.ctl")){
	die "\n\nMissing or missnamed control file: Must have the default name, i.e., \"baseml.ctl\"";
}else{
	open(CTLFILE,"baseml.ctl");
	while(<CTLFILE>){
		chomp(my $curline = $_);
		if($curline =~ /seqfile \=/){
			@seq_file_line = split(/\=/,$curline);
			$data_file_name = $seq_file_line[1];
			$data_file_name =~ s/\s+//g;
			@seq_file_line = split(/\*/,$data_file_name);
			$data_file_name = $seq_file_line[0];
			print "\n\n->datafile: $data_file_name<-\n";
		}
		elsif($curline =~ /treefile \=/){
			@tree_file_line = split(/\=/,$curline);
			$tree_file_name = $tree_file_line[1];
			$tree_file_name =~ s/\s+//g;
			@tree_file_line = split(/\*/,$tree_file_name);
			$tree_file_name = $tree_file_line[0];
			print "->treefile: $tree_file_name<-\n\n";
		}
		elsif($curline =~ / model \=/){
			@subModel_line = split(/\=/,$curline);
			$subModel = $subModel_line[1];
			$subModel =~ s/\s+//g;
			@subModel_line = split(/\*/,$subModel);
			$subModel = $subModel_line[0];
			print "->Substitution model: $subModel<-\n\n";
		}
		elsif($curline =~ /ClusMethod \=/){
			@ClusMethod_line = split(/\=/,$curline);
			$ClusMethod = $ClusMethod_line[1];
			$ClusMethod =~ s/\s+//g;
			@ClusMethod_line = split(/\*/,$ClusMethod);
			$ClusMethod = $ClusMethod_line[0];
			print "->ClusMethod: $ClusMethod<-\n\n";
		}
		elsif($curline =~ / ndata \=/){
			@ndata_line = split(/\=/,$curline);
			$ndata = $ndata_line[1];
			$ndata =~ s/\s+//g;
			@ndata_line = split(/\*/,$ndata);
			$ndata = $ndata_line[0];
			print "->ndata: $ndata<-\n\n";
		}
		elsif($curline =~ / getSE \=/){
			@ndata_line = split(/\=/,$curline);
			$get_SE = $ndata_line[1];
			$get_SE =~ s/\s+//g;
			@ndata_line = split(/\*/,$get_SE);
			$get_SE = $ndata_line[0];
			print "->get_SE: $get_SE<-\n\n";
			print "->ndata: $ndata<-\n\n";
		}
		elsif($curline =~ / sge \=/){
			@ndata_line = split(/\=/,$curline);
			$sge = $ndata_line[1];
			$sge =~ s/\s+//g;
			@ndata_line = split(/\*/,$sge);
			$sge = $ndata_line[0];
			print "->sge: $sge<-\n\n";
		}
		elsif($curline =~ / fix_kappa \=/){
			@ndata_line = split(/\=/,$curline);
			$fix_kappa = $ndata_line[1];
			$fix_kappa =~ s/\s+//g;
			@ndata_line = split(/\*/,$fix_kappa);
			$fix_kappa = $ndata_line[0];
			print "->fix_kappa: $fix_kappa<-\n\n";
		}
		elsif($curline =~ / kappa \=/){
			@ndata_line = split(/\=/,$curline);
			$kappa = $ndata_line[1];
			# $kappa =~ s/\s+//g;
			@ndata_line = split(/\*/,$kappa);
			$kappa = $ndata_line[0];
			print "->kappa: $kappa<-\n\n";
		}
		elsif($curline =~ / fix_alpha \=/){
			@ndata_line = split(/\=/,$curline);
			$fix_alpha = $ndata_line[1];
			# $fix_alpha =~ s/\s+//g;
			@ndata_line = split(/\*/,$fix_alpha);
			$fix_alpha = $ndata_line[0];
			print "->fix_alpha: $fix_alpha<-\n\n";
		}
		elsif($curline =~ / alpha \=/){
			@ndata_line = split(/\=/,$curline);
			$alpha = $ndata_line[1];
			# $alpha =~ s/\s+//g;
			@ndata_line = split(/\*/,$alpha);
			$alpha = $ndata_line[0];
			print "->alpha: $alpha<-\n\n";
		}
		elsif($curline =~ / ncatG \=/){
			@ndata_line = split(/\=/,$curline);
			$ncatG = $ndata_line[1];
			$ncatG =~ s/\s+//g;
			@ndata_line = split(/\*/,$ncatG);
			$ncatG = $ndata_line[0];
			print "->ncatG: $ncatG<-\n\n";
		}
		elsif($curline =~ / Small_Diff \=/){
			@ndata_line = split(/\=/,$curline);
			$Small_Diff = $ndata_line[1];
			$Small_Diff =~ s/\s+//g;
			@ndata_line = split(/\*/,$Small_Diff);
			$Small_Diff = $ndata_line[0];
			print "->Small_Diff: $Small_Diff<-\n\n";
		}
		elsif($curline =~ / cleandata \=/){
			@ndata_line = split(/\=/,$curline);
			$cleandata = $ndata_line[1];
			$cleandata =~ s/\s+//g;
			@ndata_line = split(/\*/,$cleandata);
			$cleandata = $ndata_line[0];
			print "->cleandata: $cleandata<-\n\n";
		}else{
			next;
		}
	}
	close(CTLFILE);
}

if($ClusMethod > -1){
	$start_at = 2;
}else{
	$start_at = 1;
}

for(my $i=$start_at; $i<$NbMaxClusters+1; $i++){
	print "\n K= $i";
	chdir "$curdir";
	mkdir("k"."$i", 0777) || warn "Cannot make initial directories: $!";
	copy("mybaseml","k$i/mybaseml") || die "Cannot copy mybaseml binary: $!";
	copy("$tree_file_name","k$i/$tree_file_name") || die "Cannot copy mybaseml tree: $!";
	copy("$data_file_name","k$i/$data_file_name") || die "Cannot copy mybaseml data: $!";
	chdir "$curdir/k$i";
	for(my $j=0; $j < $NbReps; $j++){
		if(!$SGE){
			$pm->start and next; # do the fork
		}
		chdir "$curdir/k$i";
		mkdir("rep"."$j", 0777) || warn "Cannot make rep directory $j in K=$i: $!";
		copy("mybaseml","rep$j/mybaseml") || die "Cannot copy mybaseml binary: $!";
		copy("$tree_file_name","rep$j/$tree_file_name") || die "Cannot copy mybaseml tree: $!";
		copy("$data_file_name","rep$j/$data_file_name") || die "Cannot copy mybaseml data: $!";
		chdir "$curdir/k$i/rep$j";
		write_ctl_file("$i");
		if($SGE){
			system "chmod +x mybaseml";
			system "qsub -cwd -j y -b y ./mybaseml";
		}else{
			system "chmod +x mybaseml; nice -n5 ./mybaseml";
		}
		if(!$SGE){
			$pm->finish; # do the exit in the child process
		}
	}
	unlink "$curdir/k$i/mybaseml";
}
if(!$SGE){
	$pm->wait_all_children;
}


print "\n\n\aAll done!\n";

#################################################################################
# SUBS                                                                          #
#################################################################################
sub write_ctl_file{
	my $cluster_id = shift @_;
	
	if(!$fix_kappa){
		$kappa = " ";
		for(my $i=0; $i<$ndata; $i++){
			$kappa .= random_uniform(1,1,10) . "  ";
		}
	}
	if(!$fix_alpha){
		$alpha = " ";
		for(my $i=0; $i<$ndata; $i++){
			$alpha .= random_uniform(1,.0001,10) . "  ";
		}
	}
	
	open(CTLFILE,">baseml.ctl") || die "\ncould not create ctl file for cluster $cluster_id... $!\n";
	print CTLFILE "      seqfile = $data_file_name 
     treefile = $tree_file_name

      outfile = clock.test_6_k$cluster_id.out       * main result file
        noisy = 0   * 0,1,2,3: how much rubbish on the screen
      verbose = 0   * 1: detailed output, 0: concise output
      runmode = 0   * 0: user tree;  1: semi-automatic;  2: automatic
                    * 3: StepwiseAddition; (4,5):PerturbationNNI 

        model = $subModel   * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
                    * 5:T92, 6:TN93, 7:REV, 8:UNREST, 9:REVu; 10:UNRESTu
                    
   NbClusters = $cluster_id * Number of clusters //SAB
   ClusMethod = $ClusMethod * 0: k-means, 1: k-medoids, -1: original method //SAB
          sge = $sge * 0: to run on most systems; 1: to run on my SGE queue
        
        Mgene = 0   * 0:rates, 1:separate; 2:diff pi, 3:diff kapa, 4:all diff

        ndata = $ndata
        clock = 6   * 0:no clock, 1:clock; 2:local clock; 3:CombinedAnalysis
    fix_kappa = $fix_kappa   * 0: estimate kappa; 1: fix kappa at value below
        kappa = $kappa  * initial or fixed kappa

    fix_alpha = $fix_alpha   * 0: estimate alpha; 1: fix alpha at value below
        alpha = $alpha   * initial or fixed alpha, 0:infinity (constant rate)
       Malpha = 0   * 1: different alpha's for genes, 0: one alpha
        ncatG = $ncatG   * # of categories in the dG, AdG, or nparK models of rates
        nparK = 0   * rate-class models. 1:rK, 2:rK&fK, 3:rK&MK(1/K), 4:rK&MK 

        nhomo = 0   * 0 & 1: homogeneous, 2: kappa for branches, 3: N1, 4: N2
        getSE = $get_SE   * 0: don't want them, 1: want S.E.s of estimates
 RateAncestor = 0   * (0,1,2): rates (alpha>0) or ancestral states

   Small_Diff = $Small_Diff
    cleandata = $cleandata  * remove sites with ambiguity data (1:yes, 0:no)?
*        icode = 0  * (with RateAncestor=1. try GC in data,model=4,Mgene=4)
   fix_blength = -1  * 0: ignore, -1: random, 1: initial, 2: fixed
        method = 0  * 0: simultaneous; 1: one branch at a time
";
	close(CTLFILE)
}
#################################################################################
sub get_pwd{
	my $OStype = ''; my($this_pwd,@dum);
	# test if OS is Windows
	system "ver > ostype.txt";
	open(TESTFILE,"ostype.txt") || warn "\ncould not open ostype.txt(1)... $!\n";
	while(<TESTFILE>){
		chomp(my $curline = $_);
		if($curline =~ /Windows/){
			$OStype = "Windows";
			next;
		}
	}
	close TESTFILE;

	# test if OS is linux/unix, & pretty much assumes this type if you get here
	if($OStype eq ''){
		system "uname > ostype.txt";
		open(TESTFILE,"ostype.txt") || die "\ncould not open ostype.txt(L)... $!\n";
		while(<TESTFILE>){
			chomp($OStype = $_);
		}
		close TESTFILE;
	}
	unlink "ostype.txt";

	if($OStype eq "Windows"){
		system "dir > pwd.txt";
		open(PWDFILE,"pwd.txt") || die "\ncould not open pwd.txt(W)... $!\n";
		while(<PWDFILE>){
			chomp(my $curline = $_);
			if($curline =~ /Directory of/){
				@dum = split(/ of /,$curline);
				$this_pwd = $dum[1];
				next;
			}
		}
		system "cls";
	}else{
		system "pwd > pwd.txt";
		open(PWDFILE,"pwd.txt") || die "\ncould not open pwd.txt(L)... $!\n";
		while(<PWDFILE>){
			chomp($this_pwd = $_);
		}
		system "clear";
	}
	close PWDFILE;
	unlink "pwd.txt";
	$this_pwd =~ s/\\/\//g if $OStype eq "Windows";

	return($OStype,$this_pwd);
}
#################################################################################
# THAT'S ALL FOLKS                                                              #
#################################################################################
