#!/util/bin/perl 

# --- Input: HAPMIX 3 column (LOCALANC) output
# --- Output:  integer file; one line per SNP, one integer per sample

use strict;
use Symbol;
use Switch;

unless ( scalar @ARGV == 4 )  {
  print 
    "usage:  extract_output5.pl <output-dir> <numsamples> <chrom> <pop-name>\n";
  exit(1);
}

my $DIR      = $ARGV[0];             # directory where LOCALANC files are
my $nsamples = $ARGV[1];             # number of samples in data set
my $c        = $ARGV[2];             # chromosome
my $admixpop = $ARGV[3];             # data set output name

my @filehandles = ();
for(my $i=0;$i<$nsamples;$i++)  {
  my $fh = gensym;
  my $fname = "$DIR/$admixpop.LOCALANC.$i.$c";
  open($fh,$fname) || die("Unable to open $fname\n");
  $filehandles[$i] = $fh;
}

my $fout = "$DIR/$admixpop" . ".LOCALANC.SAMPLED." . $c;
open(ANCOUT, ">$fout");

my $done = 0;
while ( !$done )  {

  for(my $i=0;$i<$nsamples;$i++)  {
    my $fh = $filehandles[$i];
    my $line = "";
    if ( $line = <$fh> )  {
      my $result = sample_columns($line);
      printf ANCOUT "%1d", $result;
    }
    else  {
      $done = 1;
      last;
    }
  }

  unless ( $done == 1 )  {
    printf ANCOUT "\n";
  }

}


for(my $i=0;$i<$nsamples;$i++)  {
  close($filehandles[$i]);
}
close(ANCOUT);





sub sample_columns()  {
  my $line = shift;
  chomp($line);
  my @C = split(/[\s]+/,$line);

  # could be four columns if physical position is included
  my $offset = 0;
  if ( scalar @C == 4 )  {
    $offset = 1;
  }

  my $r = rand();
  my $anc = 0;
  if ( $r < $C[0+$offset] )  {
    $anc = 2;
  }
  elsif ( $r < $C[1+$offset] +$C[0+$offset] )  {
    $anc = 1;
  }
  else  {
    $anc = 0;
  }

  return $anc;

}
