#!/bin/perl
# 
#  ~/Users/Sanjay/Simwalk_test2/share.pl ---
# 
#  $Id: share.pl,v 1.2 2003/02/13 15:39:44 xzhou Exp xzhou $
# 

#  Author:  Xiaojun Zhou <xzhou@request.mdacc.tmc.edu>
#  Keywords: 

#  Commentary:
#  * 
#  Code:

for ($fam=1; $fam<=$ARGV[0]; $fam++) {
    print "Working on family $fam ...\n";

    #
    if ($fam < 10) {
	$ibd_file='IBD15-22.00'.$fam;
    }
    elsif ($fam < 100) {
	$ibd_file='IBD15-22.0'.$fam;
    }
    elsif ($fam< 1000) {
	$ibd_file='IBD15-22.'.$fam;
    }
    else {
	die "999 families are the most I can handle!\n";
    }

    #
    open(IN, "<$ibd_file") || die "Unable to read $ibd_file: $!\n";

    %ids=();
    %ibds=();
    %kins=();
    @ss=();
    while (<IN>) {
	if (/^\s+([0-9]+)\s+([0-9]+)\s+/) {
	    if (!exists($ids{$1})) {
		$ids{$1}=1;
	    }
	    if (!exists($ids{$2})) {
		$ids{$2}=1;
	    }
	    $key=join('-', sort{$a<=>$b}($1, $2));
	    @pots=split /\s+/, "$'";
	    $ss[1]=$pots[2];
	    $ss[6]=$pots[3];
	    $ss[11]=$pots[4];
	    for ($i=0; $i<4; $i++) {
		$_=<IN>;
		@pots=split;
		$ss[2+$i]=$pots[0];
		$ss[7+$i]=$pots[1];
		$ss[12+$i]=$pots[2];
	    }
	    $kins{$key}=$pots[3];
	    $ibds{$key}=[@ss];
	}
    }

    close IN;

    #
    open(OUT, ">>loci.out") || die "Unable to write loci.out: $!\n";
    select OUT;
    $|=1;
    foreach $k (sort{$a<=>$b}(keys(%ids))) {
	foreach $j (sort{$a<=>$b}(keys(%ids))) {
	    if ($j <= $k) { next; }
	    $key=join('-', ($k, $j));
	    if (exists($ibds{$key})) {
		#print $k, " ", $j;
		$temp=$ibds{$key};
		@pots=@$temp;
		$f1=$pots[10]+$pots[11]+$pots[13]+$pots[14];
		$f2=$pots[9]+$pots[12];
		$pi=$f2+$f1*0.5;
		$ff=$pots[11]+$pots[9];
		$mm=$pots[10]+$pots[9];
		$fmmf=$pots[13]+$pots[14]+2*$pots[12];
		printf " %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f", $f1, $f2, $pi, $ff, $mm, $fmmf; 
		print "\n";
	    }
	    else {
		die "Family $i pair $k,$j has no ibd value!\n";
	    }
	}
    }
    select STDOUT;
    close OUT;
}