snpPlot <- function(x,range=NULL, yLabels=NULL, zlim=NULL,title=NULL,colRamp=12){

     require(colorRamps)

     if (missing(x)) {
            stop("Data matrix x missing. Stopped.")
        }

        if (!is.matrix(x)) {
            x <- as.matrix(x)
        }

        if (!is.numeric(x)) {
            stop("Data matrix x must be numeric. Stopped.")
        }



     min <- min(x)
     max <- max(x)

      if( !is.null(zlim) ){
       min <- zlim[1]
       max <- zlim[2]
    }

    if( !is.null(yLabels) ){
       yLabels <- c(yLabels)
    } else {
       yLabels <- c(1:nrow(x))
   }

    if( !is.null(range) ){
      xLabels <- c(range)
    } else {
      xLabels <- c(1:ncol(x))
    }

     if( is.null(title) ){
     title <- c()
    }

 # Red and green range from 0 to 1 while Blue ranges from 1 to 0

    ColorRamp <- switch(colRamp,
                    rgb( seq(0,1,length=256),seq(0,1,length=256),seq(1,0,length=256)),
                    rgb( seq(0,1,length=256),seq(0.5,1,length=256),seq(1,0,length=256)),
                    rgb( seq(0,1,length=256),seq(0,0.5,length=256),seq(1,0,length=256)),
                    blue2green(256),
                    green2red(256),
                    blue2yellow(256),
                    cyan2yellow(256),
                    magenta2green(256),
                    matlab.like(256),
                    matlab.like2(256),
                    blue2green2red(256),
                    rgb(t(rep(c(255,255,255),256) - col2rgb(matlab.like2(256))),maxColorValue = 255),
                    rgb(t(rep(c(255,255,255),256) - col2rgb(blue2green2red(256))),maxColorValue = 255),
                    rgb(t(abs(rep(c(0,255,0),256) - col2rgb(matlab.like2(256)))),maxColorValue = 255),
                    rgb(t(abs(rep(c(0,55,0),256) - col2rgb(blue2green2red(256)))),maxColorValue = 255),
                    rgb(t(abs(rep(c(255,0,255),256) - col2rgb(matlab.like2(256)))),maxColorValue = 255),
                    rgb(t(abs(rep(c(255,0,255),256) - col2rgb(blue2green2red(256)))),maxColorValue = 255)
            )

# ColorRamp <- rgb( seq(0,1,length=256),  # Red
#                   seq(0,1,length=256),  # Green
#                   seq(1,0,length=256))  # Blue
 ColorLevels <- seq(min, max, length=length(ColorRamp))

 # Reverse Y axis
 reverse <- nrow(x) : 1
 yLabels <- yLabels[reverse]
 x <- x[reverse,]

 # Data Map
# par(mar = c(3,5,2.5,2))
# par(mar = c(3,5,2.5,1))
# Original: par(mar=c(5.1,4.1,4.1,2.1))
 marO <- par("mar")
 par(mar = c(2.5,6,3,1))
 image(1:ncol(x), 1:nrow(x), t(x), col=ColorRamp, xlab="",
 ylab="", axes=FALSE, zlim=c(min,max))
 if( !is.null(title) ){
    title(main=title)
 }

xt <- c(1,axTicks(1))

lxt <- length(xt)

if (abs(xt[lxt]-ncol(x))>5) {
    xt <- c(xt,ncol(x))
}

axis(BELOW<-1, at=xt, labels=prettyNum(xLabels[xt], big.mark = ","), cex.axis=0.7)

axis(LEFT <-2, at=1:length(yLabels), labels=yLabels, las= HORIZONTAL<-1,cex.axis=0.7)

 par(mar = marO)

 layout(1)
}


plotRangeL <- function(Lout,ibb,physPos=NULL,colRamp=12,val=c(0.0,2.0,1.0),chrom=1,count=0,labelsNA=NULL,prange=NULL,labelsNA1 = c("model L") ) {

    if (missing(Lout)) {
        stop("Sample SNVs 'Lout' are missing. Stopped.")
    }

    if (missing(ibb)) {
        stop("tagSNVs 'ibb' are missing. Stopped.")
    }



n <- dim(Lout)[2]
n1 <- length(ibb)
doMis <- FALSE
if (n1>4) {
    n1 <- n1-1
    doMis <- TRUE
  }

if ((!is.null(prange))&&(!is.null(physPos))) {

if (length(prange)==2) {
ma <- length(ibb[[1]])

a1 <- which(physPos>=prange[1])
a2 <- which(physPos<=prange[2])

aa <- intersect(a1,a2)

if (length(aa)>1) {

  ibb <- ibb[[1]][aa]
  physPos <- physPos[aa]

}

}

}

if (is.null(physPos)) {
    physPos <- ibb*100
}


ma <- length(ibb[[1]])


range <- ibb[[1]][1]:ibb[[1]][ma]

lr <- ibb[[1]][ma]-ibb[[1]][1]+1

physRanget <- physPos[1]:physPos[ma]

lt <- length(physRanget)

ltr <- lt%/%(lr-1)

ltt <- c(1,(1:(lr-2))*ltr,lt)

physRange <- physRanget[ltt]

physLength <- round(lt/1000)

physPosM <- (physPos[ma] + physPos[1])%/%2


mat <- matrix(val[1],(n+n1+1),lr)
mat[(n+1),] <- val[2]

for (iB in 1:n1) {

  ibb_pos <- match(ibb[[iB]],range)

  mat[(n+1+iB),ibb_pos] <- val[3]
}
if (doMis) {
   ibb_pos <- match(ibb[[n1+1]],range)
   mat[(n+1+3),ibb_pos] <- 0.35*val[3]
}

ibb_pos <- match(ibb[[1]],range)

for (j in 1:n) {
    lla <- which(Lout[range,j]>=0.1)
    llb <- setdiff(1:lr,lla)
    llc <- intersect(ibb_pos,lla)
    lla <- setdiff(lla,llc)

    mat[j,lla] <- val[2]
    mat[j,llb] <- val[1]
    mat[j,llc] <- val[3]

}

if (is.null(labelsNA)) {
    labelsNA <- as.character(1:n)
}

if (length(labelsNA1)!=n1) {
  if (n1>1) {
  ade <- rep("",(n1-1))
  labelsNA1 = c("model L",ade)
} else {
  labelsNA1 = c("model L")
}
}

labelsNA <- c(labelsNA,"",labelsNA1)

if (count<=0) {
snpPlot(mat,range=physRange,yLabels=labelsNA, title=paste("chr: ",chrom,"  ||  pos: ",prettyNum(physPosM, big.mark = ","),"  ||  length: ",physLength,"kbp  ||  #SNPs: ",ma,"  ||  #Samples: ",n,sep=""),colRamp=colRamp)
} else {
snpPlot(mat,range=physRange,yLabels=labelsNA, title=paste("count: ",count,"  ||  chr: ",chrom,"  ||  pos: ",prettyNum(physPosM, big.mark = ","),"  ||  length: ",physLength,"kbp  ||  #SNPs: ",ma,"  ||  #Samples: ",n,sep=""),colRamp=colRamp)
}


}


plotAncientGenomes <- function(IBD,filename,ibdC) {


samp <- IBD[[ibdC]]$noSamples        
ibb <- IBD[[ibdC]]$noSnps 
ibb <- as.integer(sort.int(as.integer(unique(ibb))))
snpPositions <- IBD[[ibdC]]$snpPositions

labels_ALL <- IBD[[ibdC]]$labelSamples
  
Lout <- readSamplesSpfabia(X=filename,samples=samp,lowerB=0,upperB=1000.0)

if (!exists("posA")) {
load(file="/seppdata/sepp/linkage/kpgp/annotationAChr1.Rda")
posA <<- annotationAChr1[,2]
neaA <<- as.character(annotationAChr1[,24])
refA <<- as.character(annotationAChr1[,4])
altA <<- as.character(annotationAChr1[,5])
}

mA <- match(snpPositions,posA,nomatch=0)

gA <- mA[which(mA>0)]

l2 <- length(gA) 

neaAG <- neaA[gA] 
altAG <- altA[gA]
refAG <- refA[gA]


ibbA <- c()
if (l2>0) {
for (it in 1:l2)  {
if (toupper(neaAG[it]) == toupper(altAG[it]))  {
  ibbA <- c(ibbA,ibb[it])
}
}
}

 
if (!exists("posN")) {
load(file="/seppdata/sepp/linkage/kpgp/neanderHg19KPGPChr1.Rda")
posN <<- neanderHg19$V2
neaN <<- as.character(neanderHg19$ad1)
refN <<- as.character(neanderHg19$V5)
altN <<- as.character(neanderHg19$V6)
}

mN <- match(snpPositions,posN,nomatch=0)

gN <- mN[which(mN>0)]

l2 <- length(gN) 

neaNG <- neaN[gN] 
altNG <- altN[gN]
refNG <- refN[gN]

ibbN <- c()
ibbNa <- c()
if (l2>0) {
for (it in 1:l2)  {
  AAs <- strsplit(as.character(neaNG[it]),",")[[1]]
if (length(AAs)>1)  {
if (toupper(AAs[1]) == toupper(altNG[it])) {
  ibbN <- c(ibbN,ibb[it])
} 
} else {
  if (AAs[1]=="-") {
    ibbNa <- c(ibbNa,ibb[it])
  }
}
}
}

if (!exists("posD")) {
load(file="/seppdata/sepp/linkage/kpgp/denisovaHg19KPGPChr1.Rda")
posD <<- denisovaHg19$V2
neaD <<- as.character(denisovaHg19$ad1)
refD <<- as.character(denisovaHg19$V5)
altD <<- as.character(denisovaHg19$V6)
}

mD <- match(snpPositions,posD,nomatch=0)

gD <- mD[which(mD>0)]


l2 <- length(gD) 

neaDG <- neaD[gD]  
altDG <- altD[gD]
refDG <- refD[gD]


ibbD <- c()
if (l2>0) {
for (it in 1:l2)  {
  AAs <- strsplit(as.character(neaDG[it]),",")[[1]]
if (length(AAs)>1)  {
if (toupper(AAs[1]) == toupper(altDG[it])) {
  ibbD <- c(ibbD,ibb[it])
}
}
}
}

ibbL <- list(ibb,ibbA,ibbN,ibbD,ibbNa)
             
plotRangeL(Lout=Lout,ibb=ibbL,physPos=snpPositions,colRamp=12,val=c(0.0,2.0,1.0),chrom=1,count=0,labelsNA=labels_ALL,labelsNA1=c("model L","Ancestor","Neandertal","Denisova"))

}