#!/usr/bin/perl

use strict;
use warnings;
use Data::Dumper;
use File::Basename;
use Getopt::Long;
use constant MESSAGE => 1;

my ($fof);
GetOptions('fof=s' => \$fof);

usage() unless defined $fof;

my $fof_fh = fopen($fof);
my @phase_out = <$fof_fh>;
fclose($fof_fh);

chomp @phase_out;

my $all_nt = [];
my $all_ind = [];
foreach my $file (@phase_out) {
	my ($nt, $ind) = read_phase_output($file);
	push @$all_nt, @$nt;
	push @$all_ind, @$ind;
}
$all_nt = filter($all_nt);
output($all_nt, $all_ind);


# --- sub ---
sub message
{
	my @msg = @_;
	print STDERR join(" ", @msg), "\n" if MESSAGE;
}

sub fopen
{
	my ($file, $method) = @_;
	my ($fh);
	if (defined $method and $method eq 'w') {
		open $fh, ">", $file or die "$file: $!\n";
	} else {
		open $fh, $file or die "$file: $!\n";
	}
	message("$file opened");
	return $fh;
}

sub fclose
{
	my $fh = shift;
	close $fh;
	message("file handle closed");
}

sub usage
{
	print<<USAGE;

Usage : perl $0 --fof <fof.txt>

USAGE
	exit;
}

sub read_phase_output
{
	my $phase_output = shift;
	message("reading $phase_output");

	my $output;
	my $fh = fopen($phase_output);
	{
		local $/;
		$output = <$fh>;
	}
	fclose($fh);

	my ($ind_num, $loci_num, $best_pairs, @lines);
	($ind_num, $loci_num) = $output =~ /Number of Individuals: (\d+)\nNumber of Loci: (\d+)\n/;
	($best_pairs) = $output =~ /BEGIN BESTPAIRS1\n(.*)\nEND BESTPAIRS1/s;
	@lines = split /\n/, $best_pairs;

	if (@lines != $ind_num*3) {
		message("unequal individual number");
		die;
	}

	my (@ind, $nt, $n);
	$n = 0;
	foreach my $line (@lines) {
		if ($line =~ /^\d+ #(\S+)/) {
			push @ind, $1;
		} else {
			$line =~ s/[^acgtACGTU\?]//g;
			my @seq = split '', $line;
			if (@seq != $loci_num) {
				message("unequal loci number: $n, $ind[$n]");
				die;
			}
			$nt->[$n] = \@seq;
			$n++;
		}
	}

	message("reading $phase_output finished");
	return ($nt, \@ind);
}

sub transpose
{
	my $mat = shift;

	my $t_mat;
	my $N = scalar @$mat;
	my $L = scalar @{$mat->[0]};
	for (my $i=0; $i<$N; $i++) {
		for (my $j=0; $j<$L; $j++) {
			$t_mat->[$j]->[$i] = $mat->[$i]->[$j];
		}
	}
	return $t_mat;
}

sub filter
{
	my $tab = shift;

	$tab = transpose($tab);
	my $L = @$tab;
	my $N = @{$tab->[0]};
	my $filtered = [];
	LOCI: for (my $i=0; $i<$L; $i++) {
		my $allele = $tab->[$i]->[0];
		IND: for (my $j=1; $j<$N-2; $j++) {
			if ($tab->[$i]->[$j] ne $allele) {
				push @$filtered, $tab->[$i];
				next LOCI;
			}
		}
	}
	$filtered = transpose($filtered);

	return $filtered;
}

sub output
{
	my ($tab, $ind) = @_;

	my $N = @$tab;
	my $L = @{$tab->[0]};
	print "  ;1.0\n";
	map { printf("%05d ;", $_) } (1 .. $L);
	print "\n";
	map { print "10;" } (1 .. $L);
	print "\n";
	foreach my $name (@$ind) {
		print ">", join(";", $name."A", 1, "", "", "", ""), "\n";
		my $seq = shift @$tab;
		print join('', @$seq), "\n";
		print ">", join(";", $name."B", 1, "", "", "", ""), "\n";
		$seq = shift @$tab;
		print join('', @$seq), "\n";
	}
}


