#!/usr/bin/perl
#
use strict;
use warnings;
use Data::Dumper;
use Getopt::Std;
use constant DEBUG => 1;

our (%opts, $out_fh, $err_fh);

getopts('i:o:e:p:', \%opts);

if (!defined $opts{i}) {
	usage();
	exit;
} else {

	if (!defined $opts{e}) {
		$err_fh = *STDERR;
	} else {
		open $err_fh, ">", $opts{e} or die join(" ", $opts{e}, $!), "\n";
	}

	if (!defined $opts{o}) {
		$out_fh = *STDOUT;
	} else {
		open $out_fh, ">", $opts{o} or die join(" ", $opts{o}, $!), "\n";
	}

	if (defined $opts{p}) {
		open FH, "<", $opts{p} or die join(" ", $opts{p}, $!), "\n";
		my @pop_list = <FH>;
		chomp @pop_list;
		close FH;
		convert($opts{i}, @pop_list);
	} else {
		convert($opts{i});
	}

	close $out_fh;
	close $err_fh;
}

sub usage
{
	print <<USAGE;
Usage : $0 -i <snp_table> [-e] <error.txt> [-o] <output.txt> [-p] <population_list>
USAGE
	exit;
}

sub message
{
	my @msg = @_;
	print $err_fh "@msg", "\n" if DEBUG;
}

sub convert
{
	message('start reading snp table.');

	my ($snp_tab, @pop_list) = @_;
	my (
			$header,
			$line,
			$ind,
			$genotype,
			$population,
			$SNPs,
		); 
	my (
			@lines,
			@position,
			@individual,
			@fst_pos,
			@sec_pos,
			@genotypes,
		);

	open FH, $snp_tab or die "$snp_tab: $!\n";
	@lines = <FH>;
	close FH;

	map {$population->{$_} = 1} @pop_list;

	map {shift @lines} (1 .. 3); # skip the first 3 lines;
	$header = shift @lines;
	@position = sort { $a <=> $b } grep { /\d+/ } (split /\s+/, $header);
	if ($position[0] < 0) {
		my $offset = 0 - $position[0];
		map { $_ += $offset } @position;
	}
	@fst_pos = map {$_ * 2} (0 .. $#position);
	@sec_pos = map {$_ * 2 + 1} (0 .. $#position);

	while ($line = shift @lines) {
		chomp $line;
		#next unless $line =~ /^(\S+)\t((\t?[ACGT\?]\t[ACGT\?])+)$/;
		($ind, @genotypes) = split /\t/, $line;
		next if @pop_list and !defined $population->{$ind};

		my @missing = grep { /\?/ } @genotypes;
		next if scalar(@missing) / scalar(@genotypes) > 0.3;

		my @haptype1 = @genotypes[@fst_pos];
		my @haptype2 = @genotypes[@sec_pos];

		$SNPs->{$ind} = [\@haptype1, \@haptype2];
	}

	if (keys %$SNPs == 0) {
		return;
	}

	print $out_fh scalar(keys %$SNPs), "\n";
	print $out_fh scalar(@position), "\n";
	print $out_fh join(" ", 'P', @position), "\n";
	print $out_fh 'S' x scalar(@position), "\n";

	foreach $ind (sort { $a cmp $b } keys %$SNPs) {
		print $out_fh '#', $ind, "\n";
		print $out_fh join("", @{$SNPs->{$ind}->[0]}), "\n";
		print $out_fh join("", @{$SNPs->{$ind}->[1]}), "\n";
	}

	message('finish reading snp table.');
}
