use strict;
use warnings;

my ($in_trp_res_file, $in_trp_null_res_file, $output_file) = @ARGV;
if(@ARGV < 3)
{
	die "USAGE: perl summarize_results\.pl <IN_TRP_RES_FILE_AS_SET_IN_outFile> <IN_TRP_NULL_RES_FILE_AS_SET_IN_outFileNullParams> <OUTPUT_FILE>\n";
}

## parameters that I will read from the files ##
my $alt_pi1 = 'na';
my $alt_mu = 'na';
my $alt_r = 'na';
my $alt_p = 'na';
my $alt_tree_length = 'na';
my $alt_alpha = 'na';
my $alt_kappa = 'na';

my $null_log_like = 'na';
my $alt_log_like = 'na';

my $null_pi1 = 'na';
my $null_mu = 'na';
my $null_tree_length = 'na';
my $null_alpha = 'na';
my $null_kappa = 'na';
#################################################

parse_trp_result_file($in_trp_res_file);
parse_trp_null_result_file($in_trp_null_res_file);

my $chisq_reject_5percent_1df = 3.84;
my $chisq_reject_5percent_2df = 5.99;


open (my $out, ">", $output_file) or die "could not open $output_file for writing";

print $out "##################################################################\n\n";
print $out "TraitRateProp searched for the following parameters:\n\n";

print $out "--- Phenotypic \(character\) Model Parameters ---\n";
print $out "pi1: the rate of change from state 0 to state 1\n";
print $out "mu: a factor to adapt the branch lengths of the tree to the expected number of phenotypic changes per unit time\n\n";

print $out "--- Sequence Model parameters ---\n";
print $out "alpha: controls the discrete gamma distribution (4 rate categories) that allows rate variation across sequence sites\n";
print $out "kappa: the ratio between the rates of transitions and transversions (HKY85 model)\n\n";

print $out "--- Parameters connecting the trait and sequence evolutionary processes ---\n";
print $out "r: the ratio between the rates of sequence evolution in phenotypic state 1 and in phenotypic state 0\n";
print $out "p: the proportion of sequence sites in association with the phenotypic trait\n\n";

print $out "--- Additional parameters ---\n";
print $out "tree length: the total length of the tree is optimized by shrinking/expanding it up to a factor\n\n";
print $out "##################################################################\n\n";

print $out "\n\n";

print $out "##################################################################\n\n";
print $out "TraitRateProp found the following values for the alternative model:\n\n";
print $out "--- Phenotypic \(character\) Model Parameters ---\n";
print $out "pi1: $alt_pi1\n";
print $out "mu: $alt_mu\n\n";

print $out "--- Sequence Model parameters ---\n";
print $out "alpha: $alt_alpha\n";
if($alt_kappa ne 'na')
{
	print $out "kappa: $alt_kappa\n";
}
print "\n";

print $out "--- Parameters connecting the trait and sequence evolutionary processes ---\n";
print $out "r: $alt_r\n";
print $out "p: $alt_p\n\n";

print $out "--- Additional parameters ---\n";
print $out "tree length: $alt_tree_length\n\n";
print $out "log-likelihood: $alt_log_like\n\n";
print $out "##################################################################\n\n";


print $out "\n\n";

print $out "##################################################################\n\n";
print $out "TraitRateProp found the following values for the null model:\n\n";
print $out "--- Phenotypic \(character\) Model Parameters ---\n";
print $out "pi1: $null_pi1\n";
print $out "mu: $null_mu\n\n";

print $out "--- Sequence Model parameters ---\n";
print $out "alpha: $null_alpha\n";
if($null_kappa ne 'na')
{
	print $out "kappa: $null_kappa\n";
}
print "\n";

print $out "--- Parameters connecting the trait and sequence evolutionary processes ---\n";
print $out "These parameters are not optimized under the null model\n\n";


print $out "--- Additional parameters ---\n";
print $out "tree length: $null_tree_length\n\n";
print $out "log-likelihood: $null_log_like\n\n";
print $out "##################################################################\n\n";


print $out "\n\n";

print $out "##################################################################\n\n";
my $D_stat = 2*($alt_log_like - $null_log_like);
$D_stat = sprintf("%.4f",$D_stat);
print $out "Chi-Sq Hypothesis Testing\n\n";
print $out "Your D statistic is: $D_stat\n";
print $out "With significance level 0\.05 ";
if($D_stat < $chisq_reject_5percent_1df)
{
	print $out "the null hypothesis cannot be rejected\n";
	print $out "No association between the phenotypic trait and the rate of sequence evolution was observed\n\n";
}
if($D_stat > $chisq_reject_5percent_2df)
{
	print $out "the null hypothesis can be rejected\n";
	print $out "An association between the phenotypic trait and the rate of sequence evolution was observed\n\n";
}
if(($D_stat > $chisq_reject_5percent_1df) and ($D_stat < $chisq_reject_5percent_2df))
{
	print $out "the null hypothesis can be rejected with 1 degree of freedom\n";
	print $out "If you compared the TraitRate alternative model \(p=1\) to the null you can reject\n\n";
}
print $out "##################################################################\n\n";
close ($out);




sub parse_trp_result_file
{
	my($in_trp_res_file) = @_;
	open(my $in,"<",$in_trp_res_file) or die "$in_trp_res_file' for reading";
	my @all_res_lines = <$in>;
	chomp(@all_res_lines);
	close($in);
	
	#charModelParam1=0.286162	charModelParam2=11.176	relativeRate=2.98559	proportion=0.83875	treeLength=1.82467	alpha=0.927738	kappa=2.13214
	#LogLikelihood = -9115.46
	#LogLikelihood Model 0 = -9190.39
		
	foreach my $line (@all_res_lines)
	{
		if($line =~ m/charModelParam1=\s*(\S+)/)
		{
			$alt_pi1 = $1;
		}
		if($line =~ m/charModelParam2=\s*(\S+)/)
		{
			$alt_mu = $1;
		}
		if($line =~ m/relativeRate=\s*(\S+)/)
		{
			$alt_r = $1;
		}
		if($line =~ m/proportion=\s*(\S+)/)
		{
			$alt_p = $1;
		}
		if($line =~ m/treeLength=\s*(\S+)/)
		{
			$alt_tree_length = $1;
		}
		if($line =~ m/alpha=\s*(\S+)/)
		{
			$alt_alpha = $1;
		}
		if($line =~ m/kappa=\s*(\S+)/)
		{
			$alt_kappa = $1;
		}
		if($line =~ m/LogLikelihood\s*=\s*(-\S+)/)
		{
			$alt_log_like = $1;
		}
		if($line =~ m/LogLikelihood\s+Model\s+0\s*=\s*(-\S+)/)
		{
			$null_log_like = $1;
		}
	}
}



sub parse_trp_null_result_file
{
	my($in_trp_null_res_file) = @_;
	open(my $in,"<",$in_trp_null_res_file) or die "$in_trp_null_res_file' for reading";
	my @all_res_lines = <$in>;
	chomp(@all_res_lines);
	close($in);

	#charModelParam1=0.405079
	#charModelParam2=11.176
	#treeLength=1.82467
	#alpha=1.03082
	#kappa=2.13214
	
		
	foreach my $line (@all_res_lines)
	{
		if($line =~ m/charModelParam1=\s*(\S+)/)
		{
			$null_pi1 = $1;
		}
		if($line =~ m/charModelParam2=\s*(\S+)/)
		{
			$null_mu = $1;
		}
		if($line =~ m/treeLength=\s*(\S+)/)
		{
			$null_tree_length = $1;
		}
		if($line =~ m/alpha=\s*(\S+)/)
		{
			$null_alpha = $1;
		}
		if($line =~ m/kappa=\s*(\S+)/)
		{
			$null_kappa = $1;
		}
	}
}