R/gaDiscreteT1.R

In CNORfuzzy: Addon to CellNOptR: Fuzzy Logic

#
#  This file is part of the CNO software
#
#  Copyright (c) 2011-2012 - EBI - Massachusetts Institute of Technology
#
#  File author(s): CNO developers (cno-dev@ebi.ac.uk)
#
#  Distributed under the GPLv2 License.
#  See accompanying file LICENSE.txt or copy at
#      http://www.gnu.org/licenses/gpl-2.0.html
#
#  CNO website: http://www.ebi.ac.uk/saezrodriguez/software/cno
#
##############################################################################
#
# File author(s): M.K Morris <morris.melody@gmail.com>
gaDiscreteT1 <-
function(CNOlist, model, paramsList, initBstring=NULL, sizeFac=0.0001,
    NAFac=1, popSize=50, pMutation=0.5, maxTime=60, maxGens=500,
    stallGenMax=100, selPress=1.2,elitism=5, relTol=0.1, verbose=TRUE,maxSizeHashTable=1000)
{

    if ((class(CNOlist)=="CNOlist")==FALSE){
        CNOlist = CellNOptR::CNOlist(CNOlist)
    } 


    simList<-prep4simFuzzy(model=model,paramsList=paramsList, verbose=FALSE)
    indexList<-indexFinder(CNOlist=CNOlist,model=model, verbose=FALSE)
    if (is.null(initBstring)==TRUE){
        initBstring <- (sample.int(dim(paramsList$type2Funs)[1],
            (simList$numType1+simList$numType2),replace=TRUE)) - 1
    }


    #initialise
    bLength<-length(initBstring)

    # mkm initialization should be changed to have multiple discrete numbers
    Pop<-rbind(initBstring,round(matrix(runif(bLength*(popSize-1),min=0,max=dim(paramsList$type1Funs)[1]), nrow=(popSize-1),ncol=bLength)))
    bestbit<-Pop[1,]
    bestobj<-Inf
    stop<-FALSE
    obj<-rep(0,popSize)
    g<-0
    stallGen<-0

    res<-rbind(
        c(g,bestobj,toString(bestbit),stallGen,Inf,Inf,toString(bestbit),0),
        c(g,bestobj,toString(bestbit),stallGen,Inf,Inf,toString(bestbit),0))
    colnames(res)<-c("Generation","Best_score","Best_bitString","Stall_Generation",
        "Avg_Score_Gen","Best_score_Gen","Best_bit_Gen","Iter_time")
    PopTol<-rep(NA,bLength)
    PopTolScores<-NA
    nInTot = length(which(model$interMat==-1))

    #Function that produces the score for a specific bitstring
    getObj<-function(x, scoresHash=NULL){
        intString<-x

        # the hash table is used to speed up code. gain is guaranteed to be at least equal to elitism/popsize
        if (is.null(scoresHash)==FALSE){
            thisScore <- scoresHash[rownames(scoresHash) == paste(unlist(x), collapse=","),1]
             if (length(thisScore) != 0){
                 return(thisScore)
            } # otherwise let us keep going
        }


        Score = computeScoreFuzzy(CNOlist=CNOlist,model=model,simList=simList, 
            indexList=indexList, paramsList, intString=x, sizeFac=sizeFac, NAFac=NAFac)

        return(Score)

    }

    #Loop
    t0<-Sys.time()
    t<-t0

    # used by the scores hashTable.
    scoresHash <- data.frame()
    # if you do want the hastable, uncomment the following line.
    #scoresHash = NULL

    while(!stop){

        #compute the scores
        scores<-apply(Pop,1,getObj, scoresHash=scoresHash)

        # fill the hash table to speed up code
        scoresHash<-fillHashTable(scoresHash, scores, Pop, maxSizeHashTable)

        #Fitness assignment: ranking, linear
        rankP<-order(scores,decreasing=TRUE)
        Pop<-Pop[rankP,]
        scores<-scores[rankP]
        fitness<-2-selPress+(2*(selPress-1)*(c(1:popSize)-1)/(popSize-1))

        #selection:stochastic uniform sampling
        wheel1<-cumsum(fitness/sum(fitness))
        breaks<-runif(1)*1/popSize
        breaks<-c(breaks,breaks+((1:(popSize-1)))/popSize)
        sel<-rep(1,popSize)

        for(i in 1:length(breaks)){
            sel[i]<-which(wheel1>breaks[i])[1]
            }

        #intermediate generation
        Pop2<-Pop[sel,]
        PSize2<-dim(Pop2)[1]
        PSize3<-popSize-elitism

        #Recombination: uniform: each bit has a .5 proba of being inherited from each parent
        mates<-cbind(ceiling(runif(PSize3)*PSize2),ceiling(runif(PSize3)*PSize2))

        #This holds the probability, for each bit, to be inherited from parent 1 (if TRUE) or 2 (if FALSE)
        InhBit<-matrix(runif((PSize3*bLength)),nrow=PSize3,ncol=bLength)
        InhBit<-InhBit < 0.5

        #Try one point crossover
        #xover<-ceiling(runif(PSize3)*(bLength-1))
        #indices<-matrix(1:bLength,nrow=PSize3,ncol=bLength,byrow=TRUE)
        #InhBit<-matrix(rep(FALSE,PSize3*bLength),nrow=PSize3,ncol=bLength)
        #for(i in 1:PSize3){
        #    InhBit[i,]<-indices[i,]<xover[i]
        #    }
        #
        # mkm this was changed in the matlab function, but I think it will work for fuzzy as implemented here

        Pop3par1<-Pop2[mates[,1],]
        Pop3par2<-Pop2[mates[,2],]
        Pop3<-Pop3par2
        Pop3[InhBit]<-Pop3par1[InhBit]

        #Mutation
        #mkm this will need to be done differently
        MutProba<-matrix(runif((PSize3*bLength)),nrow=PSize3,ncol=bLength)
        MutProba<-(MutProba < (pMutation/bLength))
        Pop3[MutProba]<-(sample.int(dim(paramsList$type2Funs)[1],length(Pop3[MutProba]),replace = TRUE)) - 1
        #Compute stats
        t<-c(t,Sys.time())
        g<-g+1
        thisGenBest<-scores[length(scores)]
        thisGenBestBit<-Pop[length(scores),]

        if(is.na(thisGenBest)){
            thisGenBest<-min(scores, na.rm=TRUE)
            thisGenBestBit<-Pop[which(scores == thisGenBest)[1],]
        }

        if(thisGenBest < bestobj){
            bestobj<-thisGenBest
            bestbit<-thisGenBestBit
            stallGen<-0
            }else{
                stallGen<-stallGen+1
                }

        resThisGen<-c(
            g,
            bestobj,
            toString(bestbit),
            stallGen,
            (mean(scores,na.rm=TRUE)),
            thisGenBest,
            toString(thisGenBestBit),
            as.numeric((t[length(t)]-t[length(t)-1]), units="secs"))

        names(resThisGen)<-c("Generation","Best_score","Best_bitString","Stall_Generation",
            "Avg_Score_Gen","Best_score_Gen","Best_bit_Gen","Iter_time")

        if(verbose) print(resThisGen)

        res<-rbind(res,resThisGen)

        #Check stopping criteria
        Criteria<-c((stallGen > stallGenMax),
            (as.numeric((t[length(t)]-t[1]), units="secs") > maxTime),
            (g > maxGens))
        if(any(Criteria)) stop<-TRUE

        #Check for bitstrings that are within the tolerance of the best bitstring
        tolScore<-scores[length(scores)]*relTol
        TolBs<-which(scores <= scores[length(scores)]+tolScore)
        if(length(TolBs) > 0){
            PopTol<-rbind(PopTol,Pop[TolBs,])
            PopTolScores<-c(PopTolScores,scores[TolBs])
            }
        else {
            PopTol <- thisGenBestBit
            PopTolScored <- thisGenBest
        }
        if(elitism > 0){
            Pop<-rbind(Pop3,Pop[(popSize-elitism+1):popSize,])
            }else{
                Pop<-Pop3
                }
    }
    #end of the while loop

    PopTol<-as.matrix(PopTol[-1,])
    PopTolScores<-PopTolScores[-1]
    TolBs<-which(PopTolScores <= scores[length(scores)]+tolScore)
    PopTol<-as.matrix(PopTol[TolBs,])
    PopTolScores<-PopTolScores[TolBs]
    PopTolT<-cbind(PopTol,PopTolScores)
    PopTolT<-unique(PopTolT,MARGIN=1)

    if(!is.null(dim(PopTolT))){
        PopTol<-PopTolT[,1:(dim(PopTolT)[2]-1)]
        PopTolScores<-PopTolT[,dim(PopTolT)[2]]
        }else{
            PopTol<-PopTolT[1:(length(PopTolT)-1)]
            PopTolScores<-PopTolT[length(PopTolT)]
            }

    res<-res[3:dim(res)[1],]
    rownames(res)<-NULL
    return(list(bString=bestbit,
        bScore=bestobj,
        currBest=scores[length(scores)],
        results=res,
        stringsTol=PopTol,
        stringsTolScores=PopTolScores))
   }

# simple function to shift a data.frame
shift <- function(d, k) rbind( tail(d,k), head(d,-k), deparse.level = 0 )



fillHashTable <-function(scoresHash, scores, Pop, maxSizeHashTable=1000)
{
    # if not a data.frame, just return NULL
    if (is.null(scoresHash)==TRUE){ return(NULL)}

    popSize = dim(Pop)[1]
    for (i in 1:dim(Pop)[1]){
        thisScore <- scoresHash[rownames(scoresHash) == paste(unlist(Pop[i,]), collapse=","), 1]
        # if not already stored, store the score and corresponding bitstring
        if (length(thisScore) == 0){
            # compute a new score
            thisScore <- scores[i]
            # rename the row (latest one) of the newly added score
            if (dim(scoresHash)[1]<maxSizeHashTable){
                scoresHash <- rbind(scoresHash, c(thisScore, 1))
                j = dim(scoresHash)[1]
                row.names(scoresHash)[j] = paste(unlist(Pop[i,]), collapse=",")
            }
            else{
                # shift by -1 so that first element put in the queue (that
                # we get rid of) is last
                #indices = c(1:maxSizeHashTable-popSize*2)
                #scoresHash[indices, ] = scoresHash[order(scoresHash[indices,2], decreasing=TRUE), ]
                scoresHash = shift(scoresHash, -1)


                #scoresHash = shift(scoresHash, -1)
                # overwrite last element with the latest score and bitstring
                scoresHash[maxSizeHashTable,] = c(thisScore, 1)
                row.names(scoresHash)[maxSizeHashTable] =
                    paste(unlist(Pop[i,]), collapse=",")
            }
         }
         else {
             count = scoresHash[rownames(scoresHash) == paste(unlist(Pop[i,]), collapse=","), 2]
             scoresHash[rownames(scoresHash) == paste(unlist(Pop[i,]), collapse=","), 2] = count+1
         }
    }
    return(scoresHash)
}