R/PostProcessChain.R

In gilles-guillot/Geneland: Clustering of living organisms from geo-referenced genetic and/or morphometrics data

#' @title PostProcessChain
#' @description   Computes posterior probabilities of population membership for each pixel of the spatial domain. 
#' @usage
#'   PostProcessChain(coordinates,
#'                    path.mcmc,nxdom, nydom,burnin)
#'   @param coordinates Spatial coordinates of individuals. A matrix with 2
#'     columns and one line per individual.
#'   
#'   @param path.mcmc  Path to output files directory 
#'   
#'   @param nxdom Number of pixel for discretization of the spatial domain
#'     in the horizontal direction
#'   
#'   @param nydom Number of pixel for discretization of the spatial domain
#'     in the vertical direction
#'   
#'   @param burnin Number of iterations of the chain to throw away.
#'     WARNING : this argument should be given the number of stored
#'     iterations (and not the number of computed iterations which differ
#'                 if \code{thinning} !=1). If you have
#'     \code{nit}=100000 and \code{thinning}=100, then only 1000 iterations
#'     are stored. Then \code{burnin}=10 will throw away 10 stored
#'     iterations, namely 100*10 computed iterations.
#' @details Posterior probability of population membership for each pixel: 
#'     They  are
#'     written in an ascii file called \file{proba.pop.membership.txt}.
#'     Two first columns are coordinates of pixels then 
#'     one column per population (thus \code{npopmax} values are computed
#'                                for each pixel). 
#'     Images in each column of \file{proba.pop.membership.txt} are stored
#'     starting from the bottom left pixel.
#'     First line of \file{proba.pop.membership.txt}  = bottom left pixel ,
#'     second line of \file{proba.pop.membership.txt} = upward neighboor of
#'     the previous pixel, etc...
#'   
#'   
#'   Posterior probability of population membership for each
#'     individual:  They are written in a file named
#'       \file{proba.pop.membership.indiv.txt}. One line per individuals and
#'       \code{2+npopmax} columns per individual. Two first columns are
#'       spatial coordinates.
#'     
#'   
#'   Label of modal population for pixels and individuals: 
#'     They are written in files named \file{modal.pop.txt} and
#'     \file{modal.pop.indiv.txt}    respectively.
#' @export


PostProcessChain <- function(coordinates=NULL,#data
                             path.mcmc, # path to MCMC output directory
                             nxdom,nydom,# resolution
                             burnin # number of saved iterations of the chain to throw away at the beginning
                             )
  {
    ## #############################
    ## read parameters of MCMC run
    print("Reading MCMC parameter file")
    param <- as.matrix(read.table(paste(path.mcmc,"parameters.txt",sep="")))
    delta.coord <- as.numeric(param[param[,1]=="delta.coord",3])
    npopmax <- as.numeric(param[param[,1]=="npopmax",3])
    nb.nuclei.max <- as.numeric(param[param[,1]=="nb.nuclei.max",3])
    nit <- as.numeric(param[param[,1]=="nit",3])
    thinning <- as.numeric(param[param[,1]=="thinning",3])
    nindiv <- as.numeric(param[param[,1]=="nindiv",3])
    use.geno2 <- as.logical(param[param[,1]=="use.geno2",3])
    use.geno1 <- as.logical(param[param[,1]=="use.geno1",3])
    use.qtc <- as.logical(param[param[,1]=="use.qtc",3])
    use.ql <- as.logical(param[param[,1]=="use.ql",3])
    spatial <- as.logical(param[param[,1]=="spatial",3])
    ## path.mcmc ends with a /
    if(substring(path.mcmc,
                 first=nchar(path.mcmc),
                 last=nchar(path.mcmc)) != "/")
      {
        path.mcmc <- paste(path.mcmc,"/",sep="")
      }
    short.path <- substring(path.mcmc,first=1,last=nchar(path.mcmc)-1)
    if(!file_test("-d",short.path)) stop(paste("Directory ", path.mcmc,
                                          "does not exist."))
    
    if(nit/thinning <= burnin)
      {
        print(paste('nit=',nit))
        print(paste('thinning=',thinning))
        print(paste('Number of saved iterations (nit/thinning) =',nit/thinning))
        print(paste('burnin=',burnin))
        stop('*** burnin has to be smaller than the number of saved iterations ***')
      }
    if(use.geno2)
      {
        nloc.geno2 <- as.numeric(param[param[,1]=="nloc.geno2",3])
        allele.numbers.geno2 <- scan(paste(path.mcmc,"allele.numbers.geno2.txt",sep=""))
        dominance <- param[param[,1]=="dominance",3]
        filter.null.alleles <- as.logical(param[param[,1]=="filter.null.alleles",3])
      }else
    {
      nloc.geno2 <- 1
      allele.numbers.geno2 <- rep(-999,nloc.geno2)
    }
    if(use.geno1)
      {
        nloc.geno1 <- as.numeric(param[param[,1]=="nloc.geno1",3])
        allele.numbers.geno1 <- scan(paste(path.mcmc,"allele.numbers.geno1.txt",sep=""))
      }else
    {
      nloc.geno1 <- 1
      allele.numbers.geno1 <- rep(-999,nloc.geno1)
    }
    if(use.ql)
      {
        nql <- as.numeric(param[param[,1]=="nql",3])
        allele.numbers.ql <- scan(paste(path.mcmc,"number.levels.ql.txt",sep=""))
      }else
    {
      nql <- 1
      allele.numbers.ql <- rep(-999,nql)
    }      
    if(use.qtc)
      {
        nqtc <- as.numeric(param[param[,1]=="nqtc",3])
      }else
    {
      nqtc <- 1
      qtc <- matrix(nrow=nindiv,ncol=nqtc,data=-999)
    }
    nalt<- c(allele.numbers.geno2,allele.numbers.geno1,allele.numbers.ql)
    nalmax <- max(1,max(nalt))
    
    ## define dummy coordinates if coord are missing
    if(is.null(coordinates))
      {
        if(spatial)
          {
            stop('Please give spatial coordinates of individuals or set argument spatial to FALSE')
          }else
        {
          n.int <- ceiling(sqrt(nindiv))
          x <- rep(seq(from=0,to=1,length=n.int),n.int)
          y <- rep(seq(from=0,to=1,length=n.int),n.int)
          y <- as.vector(t(matrix(nrow=n.int,ncol=n.int,y,byrow=FALSE)))
          coordinates <- cbind(x,y)[1:nindiv,]
        }
      }else
    {
      if(ncol(coordinates) != 2)stop('matrix of coordinates does not have 2 columns')
    }
    coordinates <- as.matrix(coordinates)
   

    
    filenpop <- paste(path.mcmc,"populations.numbers.txt",sep="")
    filenpp <- paste(path.mcmc,"nuclei.numbers.txt",sep="")
    fileu <- paste(path.mcmc,"coord.nuclei.txt",sep="")
    filec <- paste(path.mcmc,"color.nuclei.txt",sep="")
    filef <- paste(path.mcmc,"frequencies.txt",sep="")
    fileperm <- paste(path.mcmc,"perm.txt",sep="")
    filedom <- paste(path.mcmc,"proba.pop.membership.txt",sep="")
    filedomperm <- paste(path.mcmc,"proba.pop.membership.perm.txt",sep="")
    filemeanqtc <- paste(path.mcmc,"mean.qtc.txt",sep="")
    filemeanf <- paste(path.mcmc,"mean.freq.txt",sep="")
    

    ## #############################
    ## estimate npop from MCMC run
    print("Estimating number of populations")
    npop <- scan(paste(path.mcmc,"populations.numbers.txt",sep=""))
    npop.est<- order(hist(npop[-(1:burnin)],
                          breaks=seq(.5,npopmax+.5,1),plot=FALSE)$counts,
                     decreasing=TRUE)[1]

    ## find pivot = index of iteration maximizing  posterior density
    ## (among iterations where npop==npop.est and iit > burnin)
    lpd <- scan(paste(path.mcmc,"log.posterior.density.txt",sep=""))
    subK <- (npop==npop.est) & ((1:(nit/thinning)) > burnin)
    pivot <- order(lpd[subK],decreasing=TRUE)[1]
    print(paste("Iteration with highest posterior density:",pivot))

    ## #############################
    ## dimension objects for Fortran
    ncolt <- nloc.geno2 + nloc.geno1 + nql
    dom <- matrix(nrow=nxdom*nydom,ncol=npopmax,data=0)
    domperm <- matrix(nrow=nxdom*nydom,ncol=npopmax,data=0)
    coorddom <- matrix(nrow=2,ncol=nxdom*nydom,data=-999)
    indvois <- numeric(nxdom*nydom)
    distvois <- numeric(nxdom*nydom)
    orderf <- orderftmp <- 1:npopmax
    u <- matrix(nrow=2,ncol=nb.nuclei.max,data=-999)
    c <- rep(times=nb.nuclei.max,-999)
    xlim <- ylim <- rep(-999,2)
    f <- fpiv <- fmean <- array(dim=c(npopmax,ncolt,nalmax),-999)
    fmean[1:npop.est,,] <- 0
    meanqv <- meanqvpiv <- matrix(nrow=npopmax,ncol=nqtc,-999)
    ninrub=0 # nb of saved  iterations of the chain to throw away at the end


    ## objects to be written on external text files
    nitsaved <- nit/thinning
    out.orderf <- matrix(nrow=nit/thinning,ncol=npopmax,data=-999)
    ## #############################    
    ## computes posterior probabilities of population membership
    ## for pixels of the grid
    print("Calling Fortran function postprocesschain2")
    out.res<- .Fortran("postprocesschain2",
                       PACKAGE="Geneland",
                       as.integer(nxdom),
                       as.integer(nydom),
                       as.integer(burnin),
                       as.integer(ninrub),
                       as.integer(npopmax),
                       as.integer(nb.nuclei.max),
                       as.integer(nindiv),
                       as.integer(nloc.geno2),
                       as.integer(nloc.geno1),
                       as.integer(nql),
                       as.integer(ncolt),
                       as.integer(nalt),
                       as.integer(nalmax),
                       as.double(xlim),
                       as.double(ylim),
                       as.double(delta.coord),
                       as.integer(nit),
                       as.integer(thinning),
                       as.character(filenpop),
                       as.character(filenpp),
                       as.character(fileu),
                       as.character(filec),
                       as.character(filef),
                       as.character(fileperm),
                       as.character(filedom),
                       as.character(filemeanqtc),
                       as.character(filemeanf),
                       as.double(t(coordinates)),
                       as.double(u),
                       as.integer(c),
                       as.double(f),
                       as.integer(pivot),
                       as.double(fpiv),
                       as.double(fmean),
                       as.double(dom),
                       as.double(coorddom),
                       as.integer(indvois),
                       as.double(distvois),
                       as.integer(orderf),
                       as.integer(orderftmp),
                       as.integer(npop.est),
                       as.integer(use.geno2),
                       as.integer(use.geno1),
                       as.integer(use.ql),
                       as.integer(use.qtc),
                       as.integer(nqtc),
                       as.double(meanqv),
                       as.double(meanqvpiv),
                       as.integer(nitsaved),
                       as.integer(out.orderf))
    print("End of Fortran function postprocesschain2")

    ## unpack and write some outputs on external text files
    out.orderf <- matrix(nrow=nit/thinning,ncol=npopmax,
                         data=out.res[[50]])
    write.table(out.orderf,file=paste(path.mcmc,"perm.txt",sep=""),
                row.names=FALSE,col.names=FALSE,append=FALSE)
    fmean <- array(dim=c(npopmax,ncolt,nalmax),data=out.res[[34]])
    for(iloc in 1:ncolt)
      {
        append <- ifelse(iloc==1,FALSE,TRUE)
        write.table(t(fmean[,iloc,]),
              paste(path.mcmc,"mean.freq.txt",sep=""),
              row.names=FALSE,col.names=FALSE,append=append)
      }
    dom <- matrix(nrow=nxdom*nydom,ncol=npopmax,data=out.res[[35]])
    coorddom <- matrix(nrow=2,ncol=nxdom*nydom,data=out.res[[36]])
    write.table(cbind(t(coorddom),dom),
                paste(path.mcmc,"proba.pop.membership.txt",sep=""),
                append=FALSE)
                
    ## end unpack and write

    coord.grid <- read.table(filedom)[,1:2]
    pmbr <- as.matrix(read.table(filedom)[,-(1:2)])
    pmp <- rep(NA,dim(pmbr)[1])
    for(ipix in 1:(dim(pmbr)[1]))
      { pmp[ipix] <-  order(pmbr[ipix,],decreasing=TRUE)[1] }

    write.table(cbind(coord.grid,pmp),
                file=paste(path.mcmc,"modal.pop.txt",sep=""),
                quote=FALSE,row.names=FALSE,col.names=FALSE)

    ## #############################
    ## prepare arrays for pppmindiv
    indvois <- numeric(nindiv)
    distvois <- numeric(nindiv)
    u <- matrix(nrow=2,ncol=nb.nuclei.max,data=-999)
    c <- rep(times=nb.nuclei.max,-999)
    pmp <- matrix(nrow=nindiv,ncol=npopmax,data=0)
    
    ## computes posterior probabilities of population membership
    ## for individuals
    out.res<- .Fortran("pppmindiv2",
                       PACKAGE="Geneland",
                       as.integer(nindiv),
                       as.double(t(coordinates)),
                       as.integer(npopmax),
                       as.integer(nb.nuclei.max),
                       as.integer(indvois),
                       as.double(distvois),
                       as.double(u),
                       as.integer(c),
                       as.double(pmp),
                       as.character(filenpop),
                       as.character(filenpp),
                       as.character(fileu),
                       as.character(filec),
                       as.character(fileperm),
                       as.integer(nit),
                       as.integer(thinning),
                       as.integer(burnin),
                       as.integer(orderf),
                       as.integer(npop.est),
                       as.integer(pivot))
    pmp <- matrix(nrow=nindiv,ncol=npopmax,data=out.res[[9]])
    mod.pop.indiv <- t(apply(pmp,1,order))[,npopmax]

    
    ## #############################
    ## write result for pixels (with coordinates)
    write.table(cbind(coordinates,pmp),
                file=paste(path.mcmc,"proba.pop.membership.indiv.txt",sep=""),
                quote=FALSE,row.names=FALSE,col.names=FALSE)
    ## result for indiv (with coordinates)
    write.table(cbind(coordinates,mod.pop.indiv),
                file=paste(path.mcmc,"modal.pop.indiv.txt",sep=""),
                quote=FALSE,row.names=FALSE,col.names=FALSE)

    ## write  post-processing  parameters 
    param <- c(paste("nxdom :",nxdom),
               paste("nydom :",nydom),
               paste("burnin :",burnin))
    write.table(param,file=paste(path.mcmc,"postprocess.parameters.txt",sep=""),
                quote=FALSE,row.names=FALSE,col.names=FALSE)
  }