R/GABi.R

In GABi: Framework for Generalized Subspace Pattern Mining

if(getRversion() >= "2.15.1")  utils::globalVariables(c("x","fitnessArgs","diffRows"))

GABi <- function(x,nSols=0,convergenceGens=40,popsize=256,mfreq=1,xfreq=0.5,maxNgens=200,keepBest=FALSE,identityThreshold=0.75,nsubpops=4,experiod=10,diffThreshold=0.9,verbose=FALSE,maxLoop=1,fitnessArgs=list(consistency=0.8,featureWeights=rowMeans(x,na.rm=TRUE)),fitnessFun=getFitnesses.entropy,featureSelFun=featureSelection.basic){

        fitnessArgs$tabu <- list(features=0,samples=0)

        row.mins <- apply(x,MARGIN=1,min,na.rm=TRUE)
        row.maxes <- apply(x,MARGIN=1,max,na.rm=TRUE)

        invariant.rows <- which(mapply(function(x,y)x==y,x=row.mins,y=row.maxes))

        diffRows <- setdiff(1:nrow(x),invariant.rows)

        low.rows <- diffRows[which(((rowMeans(x[diffRows,],na.rm=TRUE)-row.mins[diffRows])/(row.maxes[diffRows]-row.mins[diffRows]))<(1- diffThreshold))]
        hi.rows <- diffRows[which(((rowMeans(x[diffRows,],na.rm=TRUE)-row.mins[diffRows])/(row.maxes[diffRows]-row.mins[diffRows]))> diffThreshold)]

        cat(paste(length(low.rows),"constantly low features and",length(hi.rows),"constantly high\n"))
        diffRows <- setdiff(1:nrow(x),c(low.rows,hi.rows,invariant.rows))
        cat(paste(length(diffRows),"sufficiently variable features of",nrow(x),"\n"))
        if(length(diffRows)<2){
                warning("no sufficiently variable features in dataset")
                return(list())
        }

        subpopsize <- ceiling( popsize/ nsubpops)
        demes <- list(rep(NA, nsubpops))

        # define feature selection and fitness function in current environment (so internal variables passed appropriately)
        featureSelFun <- featureSelFun
        fitnessFun <- fitnessFun
        environment(featureSelFun) <- environment()
        environment(fitnessFun) <- environment()

        allSols <- list()
        loopCount <- 0

        while(!length(allSols)>nSols & loopCount<maxLoop){

        cat(paste("initialising population for GA loop",loopCount+1,"\n"))
        for(j in 1:nsubpops){
                demes[[j]] <- array(0,dim=c(subpopsize,ncol(x)))
                for(i in 1:subpopsize){
                        demes[[j]][i,sample(ncol(x),2)]<-1
                }
        }

        bestSols <- list()

        genNo <- 0
        keptSols <- list()
        for(j in 1:nsubpops){
                bestSols[[j]] <- rep(0,ncol(demes[[j]]))
        }
        convergenceCounter <- 0
        notConverged <- TRUE
        transfergen <- experiod

        fitnesshash <- hash()

        cat("evolving GA solutions\n")

        while(notConverged & (genNo<maxNgens)){
                if(verbose) cat(paste("trying generation no.",genNo,"\n"))
                genNo <- genNo+1
                all_fitnesses <- list()
                stagnancy <- rep(FALSE,nsubpops)
                for(j in 1:nsubpops){
                        if(verbose) cat("looking up previously-calculated fitnesses\n")
                        fitnesses <- rep(0,nrow(demes[[j]]))
                        chrkeys <- apply(demes[[j]],MARGIN=1,function(x)paste(as.character(x),collapse=""))
                        hfitnesses <- lapply(chrkeys,function(x)fitnesshash[[x]])
                        prevSeen <- which(!unlist(lapply(hfitnesses,is.null)))
                        if(verbose) cat(paste(length(prevSeen),"solutions previously evaluated\n"))
                        fitnesses[prevSeen] <- unlist(hfitnesses)
                        newSols <- setdiff(1:nrow(demes[[j]]),prevSeen)
                        if(length(newSols)>1){
                                if(verbose) cat(paste("evaluating new fitnesses for",length(newSols),"sols\n"))
                                fitnesses[newSols] <- apply(demes[[j]][newSols,],MARGIN=1,fitnessFun)
                                fitnesses[is.infinite(fitnesses)] <- 0
                                fitnesses[is.na(fitnesses)] <- 0
                                for(k in 1:length(newSols)){
                                        fitnesshash[[chrkeys[newSols[k]]]] <- fitnesses[[newSols[k]]]
                                }
                        }
                        if(length(newSols)==1){
                                if(verbose) cat(paste("evaluating new fitnesses for 1 new sol\n"))
                                fitnesses[newSols] <- fitnessFun(demes[[j]][newSols,])
                                fitnesses[is.infinite(fitnesses)] <- 0
                                fitnesses[is.na(fitnesses)] <- 0
                                fitnesshash[[chrkeys[newSols]]] <- fitnesses[[newSols]]
                        }
                        all_fitnesses[[j]] <- fitnesses
                        if(verbose) cat("fitnesses evaluated\n")
                        new_bestSol <- demes[[j]][which(fitnesses==max(fitnesses))[1],]
                        identical <- as.numeric(unlist(lapply(bestSols,function(x){identical(x,new_bestSol)})))
                        if(verbose) {
                                cat("similarity of bestsol to previous bestsols:")
                                cat(identical)
                                cat("\n")
                        }
                        if(sum(identical)>0){
                                stagnancy[j] <- TRUE
                        }
                        else{
                                convergenceCounter <- 0
                        }
                        bestSols[[j]] <- new_bestSol
                        if(convergenceCounter < convergenceGens){
                                retrycount <- 0
                                while(sum(fitnesses)==0){
                                        retrycount <- retrycount+1
                                        if(verbose) cat("re-initialising population\n")
                                        demes[[j]] <- array(0,dim=c(subpopsize,ncol(x)))
                                        for(i in 1:subpopsize){
                                                demes[[j]][i,sample(ncol(x),2)]==1
                                        }
                                        if(verbose) cat("looking up previously-calculated fitnesses\n")
                                        fitnesses <- rep(0,nrow(demes[[j]]))
                                        chrkeys <- apply(demes[[j]],MARGIN=1,function(x)paste(as.character(x),collapse=""))
                                        hfitnesses <- lapply(chrkeys,function(x)fitnesshash[[x]])
                                        prevSeen <- which(!unlist(lapply(hfitnesses,is.null)))
                                        if(verbose) cat(paste(length(prevSeen),"solutions previously evaluated\n"))
                                        fitnesses[prevSeen] <- unlist(hfitnesses)
                                        newSols <- setdiff(1:nrow(demes[[j]]),prevSeen)
                                        if(length(newSols)>1){
                                                if(verbose) cat(paste("evaluating new fitnesses for",length(newSols),"sols\n"))
                                                fitnesses[newSols] <- apply(demes[[j]][newSols,],MARGIN=1,fitnessFun)
                                                fitnesses[is.infinite(fitnesses)] <- 0
                                                fitnesses[is.na(fitnesses)] <- 0
                                                for(k in 1:length(newSols)){
                                                        fitnesshash[[chrkeys[newSols[k]]]] <- fitnesses[[newSols[k]]]
                                                }
                                        }
                                        if(length(newSols)==1){
                                                if(verbose) cat(paste("evaluating new fitnesses for 1 new sol\n"))
                                                fitnesses[newSols] <- fitnessFun(demes[[j]][newSols,])
                                                fitnesses[is.infinite(fitnesses)] <- 0
                                                fitnesses[is.na(fitnesses)] <- 0
                                                fitnesshash[[chrkeys[newSols]]] <- fitnesses[[newSols]]
                                        }
                                        all_fitnesses[[j]] <- fitnesses
                                        if(retrycount > 100){
                                                diffThreshold <- diffThreshold+0.1
                                                if(diffThreshold>1){
                                                        warning("cannot find biclusters")
                                                        return(list())
                                                }
                                                if(verbose) cat(paste("no biclusters, lowering diffThreshold to",diffThreshold,"\n"))
                                                diffRows <- setdiff(1:nrow(x),invariant.rows)

                                                low.rows <- diffRows[which(((rowMeans(x[diffRows,],na.rm=TRUE)-row.mins[diffRows])/(row.maxes[diffRows]-row.mins[diffRows]))<(1-diffThreshold))]
                                                hi.rows <- diffRows[which(((rowMeans(x[diffRows,],na.rm=TRUE)-row.mins[diffRows])/(row.maxes[diffRows]-row.mins[diffRows]))>diffThreshold)]
                                                if(verbose) cat(paste(length(low.rows),"constantly low features and",length(hi.rows),"constantly high\n"))
                                                diffRows <- setdiff(1:nrow(x),c(low.rows,hi.rows,invariant.rows))
                                                if(verbose) cat(paste(length(diffRows),"differentially-expressed features of",nrow(x),"\n"))
                                                retrycount <- 0
                                        }
                                }
                                if(verbose){
                                        cat("fitnesses:\n")
                                        cat(fitnesses)
                                        cat("\n")
                                }
                                f <- fitnesses/mean(fitnesses[fitnesses>0])
                                if(keepBest & sum(fitnesses)>0){
                                        if(verbose) cat("saving best solutions\n")
                                        goodsols <- which(fitnesses>0)
                                        for(i in 1:length(goodsols)){
                                                identity <- lapply(keptSols,function(x,comp){sum(colMeans(rbind(x$chr,comp))>0.5)/min(c(sum(x$chr),sum(comp)))},comp=demes[[j]][goodsols[i],])
                                                toRemove <- which(identity>identityThreshold)
                                                newSols <- list()
                                                if(length(toRemove)>0){
                                                        newSols <- keptSols[toRemove]
                                                        keptSols <- keptSols[-toRemove]
                                                }
                                                newSols[[length(newSols)+1]] <- list(fitness=fitnesses[goodsols[i]],chr=demes[[j]][goodsols[i],])
                                                removeScores <- unlist(lapply(newSols,function(x){x$fitness}))
                                                ReplaceSol <- which(removeScores==max(removeScores))[1]
                                                keptSols[[length(keptSols)+1]] <- newSols[[ReplaceSol]]
                                        }
                                }
                                if(verbose) cat("performing selection\n")
                                intpop <- fps(population=demes[[j]],fitnesses=f,elitism=TRUE)
                                if(verbose) cat("performing reproduction\n")
                                demes[[j]] <- reproduction(intpop,xfreq,mfreq/ncol(intpop),xoverpoints=1,pinvert=0,elitism=TRUE)
                        }
                        else{
                                cat(paste("population convergence after",genNo,"generations\n"))
                                notConverged <- FALSE
                        }
                }
                if(sum(as.numeric(stagnancy))==nsubpops){
                        convergenceCounter <- convergenceCounter + 1
                }
                if(genNo==transfergen) {
                        demes <- exchangeSols(demes,all_fitnesses,TRUE,FALSE)
                        transfergen <- transfergen + experiod
                }
        }

        population <- demes[[1]]
        if(nsubpops>1){
                for(j in 2:nsubpops){
                        population <- rbind(population,demes[[j]])
                }
        }
        cat("returning final output\n")
        fitnesses <- apply(population,MARGIN=1,fitnessFun)
        fitnesses[is.infinite(fitnesses)] <- 0
        fitnesses[is.na(fitnesses)] <- 0
        goodsols <- which(fitnesses>0)
        if(length(goodsols)>0){
                if(!keepBest){
                        keptSols <- list()
                        keptSols[[1]] <- list(chr=population[goodsols[1],],fitness=fitnesses[1])
                }
                for(i in 1:length(goodsols)){
                        identity <- lapply(keptSols,function(x,comp){sum(colMeans(rbind(x$chr,comp))>0.5)/min(c(sum(x$chr),sum(comp)))},comp=population[goodsols[i],])
                        toRemove <- which(identity>identityThreshold)
                        newSols <- list()
                        if(length(toRemove)>0){
                                newSols <- keptSols[toRemove]
                                keptSols <- keptSols[-toRemove]
                        }
                        newSols[[length(newSols)+1]] <- list(fitness=fitnesses[goodsols[i]],chr=population[goodsols[i],])
                        removeScores <- unlist(lapply(newSols,function(x){x$fitness}))
                        ReplaceSol <- which(removeScores==max(removeScores))[1]
                        keptSols[[length(keptSols)+1]] <- newSols[[ReplaceSol]]
                }
        }
        if(length(keptSols)>0){
                keptSols <- lapply(keptSols,function(x){list(samples=which(x$chr==1),features=featureSelFun(which(x$chr==1)),score=fitnessFun(x$chr))})
        }
        allSols <- c(allSols,keptSols)
        cat(paste(length(allSols),"biclusters discovered\n"))

        newTabuFeatures <- unlist(lapply(keptSols,function(x)x$features))
        newTabuSamples <- unlist(lapply(keptSols,function(x)x$samples))

        fullTabuFeatures <- union(fitnessArgs$tabu$features,newTabuFeatures)
        fullTabuSamples <- union(fitnessArgs$tabu$samples,newTabuSamples)

        fitnessArgs$tabu$features <- fullTabuFeatures
        fitnessArgs$tabu$samples <- fullTabuSamples
        cat("updating arguments to re-run biclustering loop\n")
        environment(featureSelFun) <- environment()
        environment(fitnessFun) <- environment()

        loopCount <- loopCount + 1

        }

        allSols

}