use strict;
use warnings;

if (@ARGV != 3)
{
	die "USAGE: perl print_csv_per_pos\.pl <IN_TRP_LL_PER_POS_FILE_AS_SET_IN_LLPerPositionFile> <IN_TRP_RES_FILE_AS_SET_IN_outFile> <OUTPUT_FILE>\n";
}

my ($in_LL_per_pos_file, $in_trp_res_file, $out_csv_LL_per_pos_with_BF_and_posterior) = @ARGV;

my $estimated_proportion = parse_trp_result_file($in_trp_res_file);

open (my $in, "<", $in_LL_per_pos_file) or die "could not open $in_LL_per_pos_file for reading";
my @all_lines = <$in>;
chomp (@all_lines);
close ($in);

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

my $header_old = shift (@all_lines);
my @header_parts = split (/\t/, $header_old);
push (@header_parts, "empirical_posterior_Bayes");
my $header_new = join (",", @header_parts);
print $out "$header_new\n";

foreach my $line (@all_lines)
{
	# pos_index	likelihood_no_rate	likelihood_rate	bayes_factor
	my @line_parts = split (/\t/, $line);
	
	my $like_rate_not_assoc = $line_parts[1];
	my $like_rate_assoc = $line_parts[2];
	
	my $emp_post_bayes = ($like_rate_assoc * $estimated_proportion) / (($like_rate_assoc * $estimated_proportion) + ($like_rate_not_assoc * (1 - $estimated_proportion)));
	
	push (@line_parts, $emp_post_bayes);
	my $line_new = join (",", @line_parts);
	
	print $out "$line_new\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
	
	my ($alt_pi1, $alt_mu, $alt_r, $alt_p, $alt_tree_length, $alt_alpha, $alt_kappa, $alt_log_like, $null_log_like);
	
	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;
		}
	}
	
	return ($alt_p);
}