user:algorithm:ego_r [Promethee]

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

user:algorithm:ego_r [2015/10/06 07:58] (current)
Line 1: Line 1:
 +<code r>
 +#help=Efficient Global Optimization (EGO) algorithm
 +#type=Optimization
 +#output=Optimum
 +#options=initBatchSize=10,batchSize=10,iterations=10,nugget.estim=false,nugget=0.0,trend=y~1,covtype=matern5_2,liar=min,Xbounds=false,search_min=true
 +#require=lhs,DiceKriging,DiceOptim,DiceView
  
 +# iteration index
 +iEGO <<- 0
 +
 +#' constructor and initializer of R session
 +#' @returnType NULL
 +#' @return null
 +#' @author richet
 +init <- function() {
 + library(lhs)
 + library(MASS)
 + library(rgenoud)
 + library(DiceKriging)
 + library(DiceOptim)
 + library(DiceView)
 +
 +# all parameters are initialy strings, so you have to put as global non-string values
 + initBatchSize <<- as.integer(initBatchSize)
 + batchSize <<- as.integer(batchSize)
 + iterations <<- as.integer(iterations)
 + #covtype
 + nugget.estim <<- as.logical(nugget.estim)
 + nugget <<- as.numeric(nugget)
 + if (nugget == 0) {
 + nugget <<- NULL
 + }
 + trend <<- as.formula(trend)
 + if (liar == "min") {
 + liar <<- as.function(min)
 + }
 + if (liar == "max") {
 + liar <<- as.function(max)
 + }
 + noise.var <<- NULL
 + search_min <<- as.logical(search_min)
 + Xbounds <<- as.logical(Xbounds)
 + iEGO <<- 0
 +}
 +
 +#' first design building. All variables are set in [0,1]. d is the dimension, or number of variables
 +#' @param d number of variables
 +#' @returnType matrix
 +#' @return next design of experiments
 +#' @author richet
 +initDesign <- function(d) {
 + set.seed(1)
 + lhs <- maximinLHS(n=initBatchSize,k=d)
 + if (Xbounds) {
 + e=c(0,1)
 + id=1
 + while(id<d){
 + e=rbind(cbind(e,0),cbind(e,1))
 + id=id+1
 + }
 + Xinit=rbind(as.matrix(lhs),as.matrix(e))
 + } else {
 + Xinit=as.matrix(lhs)
 + }
 + return(Xinit)
 +}
 +
 +#' iterated design building.
 +#' @param X data frame of current doe variables (in [0,1])
 +#' @param Y data frame of current results
 +#' @returnType data frame or matrix
 +#' @return  next doe step
 +#' @author richet
 +nextDesign <- function(X,Y) {
 + if (iEGO > iterations) return();
 +
 + d = dim(X)[2]
 + if (dim(Y)[2] == 2) {
 + noise.var <<- as.array(Y[,2])
 + } else {
 + noise.var <<- NULL
 + }
 +
 + if (search_min) {y=Y[,1]} else {y=-Y[,1]}
 + kmi <- km(control=list(trace=FALSE),trend,optim.method='gen',penalty = NULL,covtype=covtype,nugget.estim = nugget.estim, nugget = nugget, noise.var = noise.var,design=X,response=y)
 +
 + EGOi <- max_qEI.CL(model=kmi,npoints=batchSize,L=liar(as.array(Y[,1])),lower=rep(0,d),upper=rep(1,d),control=list(trace=FALSE))
 +
 + Xnext <- EGOi$par
 +
 + iEGO <<- iEGO + 1
 + return(as.matrix(Xnext))
 +}
 +
 +#' final analysis. All variables are set in [0,1]. Return HTML string
 +#' @param X data frame of doe variables (in [0,1])
 +#'q @param Y data frame of results
 +#' @returnType String
 +#' @return HTML string of analysis
 +#' @author richet
 +analyseDesign <- function(X,Y) {
 + analyse.files <<- paste("sectionview_",iEGO-1,".png",sep="")
 + resolution <- 600
 + if (dim(Y)[2] == 2) {
 + noise.var <<- as.array(Y[,2])
 + } else {
 + noise.var <<- NULL
 + }
 + if (search_min) {
 + m = min(Y)
 + x = as.matrix(X)[Y==m,]
 + html=paste(sep="<br/>",paste("<HTML>minimum is ",m),paste(sep="","found at ",paste(collapse="= ",capture.output(x)),"<br/><img src='",analyse.files,"' width='",resolution,"' height='",resolution,"'/></HTML>"))
 + plot=paste("<Plot1D name='min'>",m,"</Plot1D>")
 +
 + d = dim(X)[2]
 + if (d == 1) {
 + plotx=paste("<Plot1D name='argmin'>",paste(x),"</Plot1D>")
 + } else if (d == 2) {
 + plotx=paste("<Plot2D name='argmin'>(",paste(collapse=",",x),")</Plot2D>")
 + } else {
 + plotx=paste("<PlotnD name='argmin'>(",paste(collapse=",",x),")</PlotnD>")
 + }
 + } else {
 + m = max(Y)
 + x = as.matrix(X)[Y==m,]
 + html=paste(sep="<br/>",paste("<HTML>maximum is ",m),paste(sep="","found at ",paste(collapse="= ",capture.output(x)),"<br/><img src='",analyse.files,"' width='",resolution,"' height='",resolution,"'/></HTML>"))
 + plot=paste("<Plot1D name='max'>",m,"</Plot1D>")
 +
 + d = dim(X)[2]
 + if (d == 1) {
 + plotx=paste("<Plot1D name='argmax'>",paste(x),"</Plot1D>")
 + } else if (d == 2) {
 + plotx=paste("<Plot2D name='argmax'>(",paste(collapse=",",x),")</Plot2D>")
 + } else {
 + plotx=paste("<PlotnD name='argmax'>(",paste(collapse=",",x),")</PlotnD>")
 + }
 + }
 +
 + model <- km(control=list(trace=FALSE),trend,optim.method='gen',penalty = NULL,covtype=covtype,nugget.estim = nugget.estim, nugget = nugget, noise.var = noise.var,design=X,response=Y)
 +
 + png(file=analyse.files,bg="transparent",height=resolution,width = resolution)
 + sectionview(model=model,center=x)
 + dev.off()
 +
 + return(paste(html,plot,plotx))
 +}
 +
 +#' temporary analysis. All variables are set in [0,1]. Return HTML string
 +#' @param X data frame of doe variables (in [0,1])
 +#' @param Y data frame of  results
 +#' @returnType String
 +#' @return HTML string of analysis
 +#' @author richet
 +analyseDesignTmp <- function(X,Y) {
 + analyseDesign(X,Y)
 +}
 +</code>
© IRSN - All right reserved - Legal information