R/GenAlgForSubsetSelectionNoTest.R

In STPGA: Selection of Training Populations by Genetic Algorithm

GenAlgForSubsetSelectionNoTest<-function (P, ntoselect, npop=100, nelite=5, keepbest=TRUE,tabu=TRUE, tabumemsize=1, mutprob=.8, mutintensity=1, niterations=500, minitbefstop=200,niterreg=5, lambda=1e-6, plotiters = FALSE,plottype=1, errorstat="PEVMEAN2",C=NULL, mc.cores=1,  InitPop=NULL, tolconv=1e-7, Vg=NULL, Ve=NULL, Fedorov=FALSE){
  if ((ncol(P)+1)>ntoselect){warning("The algorithm does not work well with p>ntrain, perhaps use a unsupervised dimension reduction on P.")}
  
  ridgereg<-function(y=CurrentPopFuncValues ,x=designforCurrentPop,lambda=lambda){
		n<-nrow(x)
		p<-ncol(x)
		mindim<-min(p,n)
		rownames(x)<-NULL
		svdX<-svd(x, nu=mindim,nv=mindim)
		insvdnonzero<-1:mindim
		diagvecforinv<-(svdX$d[insvdnonzero])/((svdX$d[insvdnonzero])^2+lambda)
		
		coef<-tcrossprod(svdX$v%*%diag(diagvecforinv),svdX$v)%*%t(x)%*%(y-mean(y))
		
		return(list(coef=coef))
  }
  
  getsetsfedorov<-function(dummyx){
    Candidates<-rownames(P)
    npc<-min(c(ntoselect-10,ncol(P)))
    randint<-sample(1:10, 1)
if (randint==1){
    optsol<-AlgDesign::optFederov(data=P[,1:npc],nTrials=ntoselect, criterion="D")
    outf<-Candidates[optsol$rows]

} else if (randint==2){
    optsol<-AlgDesign::optFederov(data=P[,1:npc],nTrials=ntoselect, criterion="A")
    outf<-Candidates[optsol$rows]

} else if (randint==3){
    optsol<-AlgDesign::optFederov(data=P[,1:npc],nTrials=ntoselect, criterion="I")
    outf<-Candidates[optsol$rows]
} else {
    outf<-sample(Candidates, ntoselect)
}
    return(outf)
    
  }
  if (is.null(InitPop)){
  if (sum(errorstat %in%c("AOPT", "DOPT", "CDMEAN", "PEVMEAN", "PEVMEAN2"))>0){
   if (Fedorov) {InitPop<-lapply(1:npop, function(q){
      return(getsetsfedorov(q))
    })
   } else {
     InitPop<-lapply(1:npop, function(q){
       return(sample(rownames(P),ntoselect))
     })
   }
    
  }
  }
  
  mc.cores <- switch( Sys.info()[['sysname']],
                      Windows = {1}, 
                      Linux   = {mc.cores},
                      Darwin  = {mc.cores})
nelite=nelite-1
 if(tabu){tabumemsize=tabumemsize}else{tabumemsize=0}
 if (tabumemsize>0){
 memoryfortabu<-vector(mode="list", length=tabumemsize)
 } else {memoryfortabu=NULL}
 if (errorstat%in%c("CDMEANMM","PEVMEANMM","GAUSSMEANMM", "userMM")){
   K=solve(P)
  }
  
  linenames<-rownames(P)
  
  completeton<-function(x){if (length(x)<ntoselect){
    x<-c(x, sample(setdiff(unique(unlist(InitPop)), x), ntoselect-length(x)))
  }
    return(x)
  }
  
  if (is.null(InitPop)){InitPop<-lapply(1:npop, function(x){return(sample(linenames, ntoselect))})
  } else {
    InitPop<-lapply(InitPop, completeton)
  }
  
  if (!is.null(errorstat)){
      if (errorstat%in%c("CDMEANMM","PEVMEANMM","GAUSSMEANMM", "userMM")){
        InitPopFuncValues <- as.numeric(unlist(mclapply(InitPop, FUN = function(x) {
          do.call(errorstat, list(x,setdiff(linenames,x),P,K,lambda,C, Vg, Ve))
                  },mc.cores=mc.cores,mc.preschedule=T)))
              } else {
    InitPopFuncValues <- as.numeric(unlist(mclapply(InitPop, FUN = function(x) {
      do.call(errorstat, list(x,setdiff(linenames,x),P,lambda, C))
    },mc.cores=mc.cores,mc.preschedule=T)))
              }
    }
  designforInitPop<-Reduce('rbind',lapply(InitPop,function(x){as.numeric(linenames%in%x)}))
  #designforInitPop<-cbind(rep(1,nrow(designforInitPop)),designforInitPop)
  betahat=ridgereg(y=InitPopFuncValues ,x=designforInitPop,lambda=lambda)$coef
      
 # betahat=solve(crossprod(designforInitPop)+lambda*diag(ncol(designforInitPop)),t(designforInitPop)%*%InitPopFuncValues)
  orderofInitPop <- order(InitPopFuncValues, decreasing = FALSE)
  ElitePop <- mclapply(orderofInitPop[1:nelite], FUN = function(x) {
    return(InitPop[[x]])
  },mc.cores=mc.cores,mc.preschedule=T)
  ElitePopFuncValues <- InitPopFuncValues[orderofInitPop[1:nelite]]
  
  orderofBetasforCurrentPop <- order(betahat, decreasing = FALSE)
  bestpredicted<-linenames[orderofBetasforCurrentPop[1:ntoselect]]
  ElitePop[[nelite+1]]<-bestpredicted
  ElitePopFuncValues[nelite+1]<-mean(ElitePopFuncValues)
  minvec <- c()
  sdvec<-c()
  designforElitePop<-Reduce('rbind',lapply(ElitePop,function(x){as.numeric(linenames%in%x)}))
  if (length(ElitePop)>2){sdofelitedesigns<-mean(apply(designforElitePop,2,sd))} else {sdofelitedesigns=0}
  sdvec<-c(sdvec, sdofelitedesigns)
  minvec <- c(minvec, min(ElitePopFuncValues))

  for (iters in 1:niterations) {
    if (iters>minitbefstop){
      maxmeans<-max(minvec[(length(minvec)-minitbefstop):length(minvec)])
      minmeans<-min(minvec[(length(minvec)-minitbefstop):length(minvec)])
      meansdiff<-maxmeans-minmeans
       if(meansdiff<tolconv){break()
     }
    }
    
    if (iters>tabumemsize){
    CurrentPop <- GenerateCrossesfromElites(Elites = ElitePop, 
                                            Candidates = linenames, npop = npop, mutprob = mutprob,mc.cores=mc.cores, mutintensity=mutintensity, memoryfortabu=memoryfortabu)
    } else {CurrentPop <- GenerateCrossesfromElites(Elites = ElitePop, 
                                                         Candidates = linenames, npop = npop, mutprob = mutprob,mc.cores=mc.cores, mutintensity=mutintensity, memoryfortabu=NULL)}
    if (tabumemsize>0){
      memoryfortabu[[iters%%tabumemsize+1]]<-CurrentPop
    } else {memoryfortabu=NULL}
    if (keepbest){CurrentPop<-c(CurrentPop, ElitePop[1])}
    
      if (!is.null(errorstat)){
        if (errorstat%in%c("CDMEANMM","PEVMEANMM","GAUSSMEANMM", "userMM")){
         CurrentPopFuncValues <- as.numeric(unlist(mclapply(CurrentPop, 
                                                         FUN = function(x) {
                                                           do.call(errorstat, list(x,setdiff(linenames,x),P,K,lambda,C, Vg, Ve))
                                                                                                                    },mc.cores=mc.cores,mc.preschedule=T))) } else {
     
        CurrentPopFuncValues <- as.numeric(unlist(mclapply(CurrentPop, 
                                                             FUN = function(x) {
                                                               do.call(errorstat,list(x, setdiff(linenames,x), P, lambda, C))
                                                             },mc.cores=mc.cores,mc.preschedule=T)))
        }
      }
  
    if (niterreg>iters){
    designforCurrentPop<-Reduce('rbind',lapply(CurrentPop,function(x){as.numeric(linenames%in%x)}))
    #designforCurrentPop<-cbind(rep(1,nrow(designforCurrentPop)),designforCurrentPop)
   
   betahat=ridgereg(y=CurrentPopFuncValues ,x=designforCurrentPop,lambda=lambda)$coef

   # betahat=solve(crossprod(designforCurrentPop)+lambda*diag(ncol(designforCurrentPop)),t(designforCurrentPop)%*%CurrentPopFuncValues)
    orderofBetasforCurrentPop <- order(betahat, decreasing = FALSE)
    bestpredicted<-linenames[orderofBetasforCurrentPop[1:ntoselect]]
    
    orderofCurrentPop <- order(CurrentPopFuncValues, decreasing = FALSE)
     ElitePop <- lapply(orderofCurrentPop[1:nelite], FUN = function(x) {
      return(CurrentPop[[x]])
    })
     
    designforElitePop<-Reduce('rbind',lapply(ElitePop,function(x){as.numeric(linenames%in%x)}))
    if (length(ElitePop)>2){sdofelitedesigns<-mean(apply(designforElitePop,2,sd))} else {sdofelitedesigns=0}
    sdvec<-c(sdvec, sdofelitedesigns)
    
    ElitePopFuncValues <- CurrentPopFuncValues[orderofCurrentPop[1:nelite]]
    ElitePop[[nelite+1]]<-bestpredicted
    ElitePopFuncValues[nelite+1]<-mean(ElitePopFuncValues)
    minvec <- c(minvec, min(ElitePopFuncValues))
    } else{
      orderofCurrentPop <- order(CurrentPopFuncValues, decreasing = FALSE)
      ElitePop <- lapply(orderofCurrentPop[1:(nelite+1)], FUN = function(x) {
        return(CurrentPop[[x]])
      })
      
      designforElitePop<-Reduce('rbind',lapply(ElitePop,function(x){as.numeric(linenames%in%x)}))
      if (length(ElitePop)>2){sdofelitedesigns<-mean(apply(designforElitePop,2,sd))} else {sdofelitedesigns=0}
      sdvec<-c(sdvec, sdofelitedesigns)
      ElitePopFuncValues <- CurrentPopFuncValues[orderofCurrentPop[1:(nelite+1)]]
      minvec <- c(minvec, min(ElitePopFuncValues))
    }  
    if (plotiters) {
        if (plottype==1){
		    par(mfrow=c(1,1))
        plot(minvec, type="b")
        abline(h=min(minvec), col="red")
        mtext(side=3, line=-3, col="blue", text=paste("\n     Minimum value of objective function is ", round(min(ElitePopFuncValues), 5), ".", sep=""), adj=0, outer=T) 

      	} else if (plottype==2){
        par(mfrow=c(2,1))
        plot(minvec, type="b")
        abline(h=min(minvec), col="red")
        mtext(side=3, line=-3, col="blue", text=paste("\n     Minimum value of objective function is ", round(min(ElitePopFuncValues), 5), ".", sep=""), adj=0, outer=T) 
        barplot(round(apply(designforElitePop,2,mean),3))
        par(mfrow=c(1,1))

      	} else if (plottype==3){
      		  par(mfrow=c(3,1))
        plot(minvec, type="b")
        abline(h=min(minvec), col="red")
        mtext(side=3, line=-3, col="blue", text=paste("\n    Divergence: ", round(sdofelitedesigns,5),".\n     Minimum value of objective function is ", round(min(ElitePopFuncValues), 5), ".", sep=""), adj=0, outer=T) 
        barplot(round(apply(designforElitePop,2,mean),3))
        plot(sdvec, type="b", col=(sdvec==0)+1)
        par(mfrow=c(1,1))

      	} else {}
      
    }
  }
  ElitePop[[nelite+2]]<-minvec
  for (i in 1:(nelite+1)){
    names(ElitePop)[i]<-paste("Solution with rank ", i, sep="")
  }
  names(ElitePop)[nelite+2]<-"Best criterion values over iterarions"
  return(ElitePop)
}