#
# read summary map and extract requested markers 
#
DATA=/var/www/Mapview
AWK=gawk

if [ "$1" = "" ] || [ "$1" = "-h" ] || [ "$1" = "--help" ] 
then
        echo "Usage: mastermap (<sib-pair>.in [ordered.in]|"
        echo "                  -merlin <merlin>.dat [ordered.map]) "
        echo "                  -merlinmap <merlin>.map [ordered.map]) "
        echo "                  -mlt2 <pedtool>.mlt2 [ordered.mlt2]) "
        echo "       Extracts positions of markers described in sib-pair"
        echo "       script or pedtool mlt2 script from master map."
        echo "       Creates a reordered script, if requested."
        echo "       e.g. mastermap chrom1.in ordered1.in"
        exit
fi
TYP="sib-pair"
switch=1
while [ $switch = 1 ]
do case "$1" in
  -mlt2) TYP="mlt2"; shift ;;
  -merlin) TYP="merlin"; shift ;;
  -merlinmap) TYP="merlinmap"; shift ;;
  -*) echo "Flag not recognised: $1" ; shift ;;
  *) switch=0 ;;
  esac
done

if ! [ -f $1 ] 
then
      grep $1 $DATA/master_map.dat | $AWK '
BEGIN{print "Locus      Chr  Pos (Mb)  Ave cM"
}
     { printf "%10s %3s   %7.3f  %6s\n",  \
               $1, $7, ($19)/1000000, $20
}'
      exit
fi

$AWK -v typ=$TYP -v loci=$1 -v fil=order$$ '

BEGIN{OFS="\t"
      if (typ=="sib-pair") {
        while (getline < loci) {
          if($1=="set" && substr($2,1,3)=="loc" && \
             (substr($4,1,3)=="mar" || substr($4,1,3)=="xma")) {
            ldbloci[$3]++ 
          }
        }
      }else if (typ=="merlin") {
        while (getline < loci) {
          if($1=="M") {
            ldbloci[$2]++ 
          }
        }
      }else if (typ=="merlinmap") {
        while (getline < loci) {
          ldbloci[$2]++ 
        }
      }else if (typ=="mlt2") {
        while (getline < loci) {
          if($3!="ZYGOSITY" && $3!="QUANTITATIVE" && $3!="QUANT" && $3!="AFF") {
            ldbloci[$1]++ 
          }
        }
      }
      delete ldbloci[""]
      delete ldbloci[" "]
      nmar=0
}
$1 in ldbloci || $2 in ldbloci || $3 in ldbloci {
      delete ldbloci[$1]
      delete ldbloci[$2]
      delete ldbloci[$3]
      nmar++
      locus[nmar]=$1
      comp[nmar]=$19
      sex_av[nmar]=$20
      chr[nmar]=$7
}
END { print "Locus      Chr  Pos (Mb)  Ave cM   Theta"
      clo=0
      last_chr=""
      for(i=1;i<=nmar;i++) {
        d=0.04*(sex_av[i]-clo) 
        if(d>=0){ theta=substr(0.5*(exp(d)-1)/(exp(d)+1),1,4) }
        if (chr[i]!=last_chr) theta=0.5
        printf "%10s %3s   %7.3f  %6s   %5.2f\n",  \
               locus[i],chr[i],comp[i]/1000000,sex_av[i], theta
      
        print locus[i],chr[i],comp[i],sex_av[i] > fil 
        clo=sex_av[i]
        last_chr=chr[i]
      }
      printf "\nUnable to match:"
      n=7
      for(i in ldbloci) {
        n++
        if (n>7) { n=1; printf "\n" }
        printf " " i
      }
      printf "\n"
}' $DATA/master_map.dat

if [ "$2" != "" ]
then
        $AWK '
        BEGIN { missing="."
                nmar=0
                nmar2=0
        }
              { nmar++
                loc[nmar]=$1; chr[nmar]=$2; pos[nmar]=$3; link[nmar]=$4
        }
        END   { i=1
                slope=1
                while (i<=nmar) {
                  while(link[i]!=missing && i <=nmar) { 
                    d1=pos[i]   
                    x1=link[i]  
                    i++ 
                  }
        	  sta=i
                  while(link[i]==missing && i<=nmar) { 
                    i++ 
                  }
        	  fin=i-1
                  d2=pos[i] ; x2=link[i]
                  if ((d2-d1)>0) slope=(x2-x1)/(d2-d1) 
        	  for(i=sta; i<=fin; i++) {
        	    link[i]=x1+slope*(pos[i]-d1)
        	  }
        	}
                for(i=1;i<=nmar;i++) {
                  print loc[i], pos[i], link[i], chr[i]
                }
        }' order$$ | sort -n -k 4,4 -k 3,3 > tmp$$
        mv tmp$$ order$$

        $AWK -v typ=$TYP -v fil=order$$ '
        BEGIN { namr=0
                nmar2=0
                while(getline < fil) {
                  nmar++
                  loc[nmar]=$1; pos[nmar]=$2; link[nmar]=$3
                }
                if (typ=="merlin" || typ=="merlinmap") print "CHR MARKER POS"
        }
              { if (typ=="sib-pair") {
                  if($1=="set" && substr($2,1,3)=="loc" && \
                     substr($4,1,3)=="mar") {
                    nmar2++
                    if (link[nmar2]!="") {
                      if (nmar2>1 && pos[nmar2]<pos[nmar2-1]) {
                        print "! Consensus and genetic orders differ"
                      }
                      printf "%s %s %s %s %6s  (%6s)\n", \
                        $1,$2,loc[nmar2],$4,link[nmar2],pos[nmar2]
                    }
                  }else{
                    print $0
                  }
                }else if (typ=="merlin") {
                  if ($1=="M") {
                    nmar2++
                    if (link[nmar2]!="") {
                      printf "%s %s %s\n", \
                        "1",loc[nmar2],link[nmar2]
                    }
                  }
                }else if (typ=="merlinmap") {
                  nmar2++
                  if (link[nmar2]!="") {
                    printf "%s %s %s\n", \
                      "1",loc[nmar2],link[nmar2]
                  }
                }else if (typ=="mlt2") {
                  if($3!="ZYGOSITY" && $3!="QUANTITATIVE" && $3!="QUANT" && $3!="AFF") {
                    nmar2++
                    if (link[nmar2]!="") {
                      printf "%s %s %s %s\n", \
                        loc[nmar2],$2,link[nmar2],loc[nmar2]
                    }
                  }else{
                    print $0
                  }
                }
        }' $1 > $2

fi

rm -f order$$