#!/usr/bin/perl -w

# Copyright 2013, Naoki Takebayashi <ntakebayashi@alaska.edu>
#
# This program 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.
#
# This program 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 this program. If not, see <http://www.gnu.org/licenses/>.

# Version: 20130612

my $usage = "\nUsage: $0 [-h] nexusFile\n".
    " -h: help\n".
    "Takes one nexus file, run model-test, and print out a paup-command \n".
    "file, which incorporate the model selected based on AIC.\n";

our ($opt_h);

use File::Basename;
use Getopt::Std;
getopts('h') || die "$usage\n";
die "$usage\n" if (defined($opt_h));

die "ERROR: Please give only 1 file\n" if (@ARGV > 1); 

my $modeltestBin = "modeltest3.7";
my $mtBlock = "/usr/share/doc/modeltest-3.7/modelblockPAUPb10";

my $criterion = 'AIC'; # other choice is 'hLRT'

#@ARGV = ('-') unless @ARGV; # take STDIN when no arg
my $nex = shift;
$base = basename($nex);

if ($nex =~ /\.nex$/) {
    $base =~ s/\.nex$//;
} elsif ($nex =~ /\.nx$/) {
    $base =~ s/\.nx$//;
}

system("mkdir -p modeltest.$base");

# run modeltest
system("cd modeltest.$base; echo \"execute ../$nex; execute $mtBlock; quit WarnTSave=no;\" | paup 1>&2; $modeltestBin < model.scores > modeltest.out");

print "#NEXUS\nBEGIN PAUP;\nexecute $nex;\nset criterion=likelihood;\n";

my $encounter = ($criterion eq 'hLRT') ? 1 : 2;
my ($cnt1, $cnt2) = (0,0);

open MT, "<modeltest.$base/modeltest.out" || 
    die "Can't open the model test out\n";
while(<MT>) {
    if( /^\s*Model selected:\s*(.+)\s*/) {
	$cnt1++;
	if ($cnt1 == $encounter) {
	    print "[! With $criterion, model selected: $1 ]\n" 
	}
    }
    if (/^Lset/) {
	$cnt2++;
	if ($cnt2 == $encounter)  {
	    print "\n$_\n\n";
	}
    }
}
close(MT);

my $template = <<"EOF";
[ exhaustive search
alltrees ; ]

[ 
outgroup '';
]

[ heuristic search ]
hsearch addseq=random nreps=10 rearrlimit=40000 limitperrep=yes;

[ bootstrap nreps=50 treefile=boot.tre search=heuristic /addseq=random rearrlimit=3000 nreps=2 limitperrep=yes; ]

[ save consensus tree
contree all/strict=no majrule=yes treefile=consensus.tre;
]

savetrees file=tree.tre root=y brlens=y  format=NEXUS;
savetrees file=tree.phy root=y brlens=y  format=phylip;
[ saveassum file=paup.assm; ]

[ write brlen file ]

[
log start file=brlens;
describetrees /brlens=yes;
log stop;
]

End;
EOF

print "$template\n";

exit;