#
# Implement Goncalo's plot of ibs mean v. variance for various
# classes of relative pair
#
read.linkped <- function (fil, cols=NULL, na.strings=c("0","x")) {
   ped <- read.table(fil, head=F, na.strings=na.strings)
   if (!is.null(cols)) {
     ped <- ped[,c(1:5,cols)]
   }
   nloc<-(ncol(ped)-5)/2
   names(ped)<-c("ped","id","fa","mo","sex",
                 paste("loc",paste(as.integer(gl(nloc,2)),c("a","b"),
                        sep="."),sep=""))
   ped
}
#
# As Fortran call
#
calc.ibs <- function (ped, rel.list=1:nrow(ped)) {
   dyn.load("calcibs.so")
   nrel<-length(rel.list)
   if (nrel<2) {
     cat("ERROR: Must be at least two relatives")
     return
   }else if (nrel>nrow(ped)) {
     cat("ERROR: Too many relatives")
     return
   }
   npairs<-as.integer(nrel*(nrel-1)/2)
   ids<-paste(ped[,1],ped[,2],sep="-")
   pair.ids<-matrix(kronecker(ids,ids,paste, sep=","),nr=length(ids))
   pair.ids<-pair.ids[!lower.tri(pair.ids,diag=T)]
   gtps<-data.matrix(ped)
   gtps[is.na(gtps)]<-0
   res<-clcibs(gtps,nrel,npairs,rel.list)
   miss<-(res$mibs==0 & res$vibs==0)
   res$mibs[miss]<-NA
   res$vibs[miss]<-NA
   return(pair.ids=pair.ids, 
          relate=res$relate, mibs=res$mibs, 
          vibs=res$vibs, typed=res$typloc)
}
clcibs <- function(ped,nrel,npairs,rel.list) {
     .Fortran("clcibs",
               nobs = as.integer(length(rel.list)),
               ncol = as.integer(ncol(ped)),
               ped = data.matrix(ped[rel.list,]),
               nrel = as.integer(nrel),
               npairs = as.integer(npairs),
               relist = as.integer(rel.list),
               relate = integer(npairs),
               mibs = double(npairs),
               vibs = double(npairs),
               typloc = integer(npairs))
}
print.ibs <- function(ibs.data) {
   as.data.frame(ibs.data)
}
#
# "Abecasis" plot
#
plot.ibs <- function(ibs.data, show.rels=1:3, crit=NULL, 
                     paint.rels=2, paint.peds=NULL) {
   visible<-(ibs.data$relate %in% show.rels)
   if (sum(visible)==0) {
     cat("No relative pairs meeting your criteria\n")
     return
   }
   mean.typed<-mean(ibs.data$typed)
   v0<-mean(sqrt(ibs.data$vibs[ibs.data$relate==0]),na.rm=TRUE)
   m0<-mean(ibs.data$mibs[ibs.data$relate==0],na.rm=TRUE)
   qm<-quantile(x$mibs[ibs.data$relate==0],probs=c(0.05,0.5,0.95),na.rm=T)

   if (is.null(crit)) {
     crit<-max(0.01,1/length(x$mibs[!is.na(x$mibs) &
                             ibs.data$relate %in% paint.rels]))
   }
   rel.colours<-c("lightgoldenrod","black","red","green")
   rel.name<-  c("Unrel","Par-Off","Siblings","Half-sibs") 

   plot(sqrt(ibs.data$vibs[visible]),
        ibs.data$mibs[visible],
        xlim=c(0,0.75),ylim=c(min(ibs.data$mibs[visible],
                               na.rm=T),1),
        xlab="Standard Deviation IBS sharing",
        ylab="Mean IBS sharing",
        main=paste("Relative pair IBS sharing over ",
                    round(mean.typed,digits=1),
                    " markers", sep=""),
        pch=1,col=rel.colours[1+ibs.data$relate[visible]])
        legend(0.6, 0.9,
               rel.name[1+show.rels],
               fill=rel.colours[1+show.rels])
   points(v0,m0, pch=3, cex=3)
   abline(h=qm[c(1,3)], lty=2)
   abline(h=qm[2], lty=3)
   text(0.4,qm[2],"Median IBS sharing for unrelated individuals",c(0,0.05))
#
# Paint the outliers
#
   if (!is.null(paint.rels)) {
     paint.qtle<-quantile(x$mibs[ibs.data$relate %in% paint.rels],
                          probs=c(crit,1-crit),na.rm=T)
     visible<-(ibs.data$relate %in% paint.rels & 
              (x$mibs<paint.qtle[1] | x$mibs>paint.qtle[2] | x$vibs<0.01))
     points(sqrt(ibs.data$vibs[visible]), ibs.data$mibs[visible],
            pch=16, col=rel.colours[1+paint.rels])
     text(sqrt(ibs.data$vibs[visible]), ibs.data$mibs[visible], 
          ibs.data$pair.ids[visible], pos=4, col=rel.colours[1+paint.rels])
   }
#
# Paint members of particular families
#
   if (!is.null(paint.peds)) {
     ped<-matrix(unlist(strsplit(ibs.data$pair.ids,"[-,]")),byr=T,nc=4)
     visible<-(ped[,1]==ped[,3] & ped[,1] %in% paint.peds)
     points(sqrt(ibs.data$vibs[visible]), ibs.data$mibs[visible],
            pch=16, col=rel.colours[1+ibs.data$relate[visible]])
     text(sqrt(ibs.data$vibs[visible]), ibs.data$mibs[visible], 
          ibs.data$pair.ids[visible], pos=4,
          col=rel.colours[1+ibs.data$relate[visible]])
   }

}