R/clusterKriging.R

In CEGO: Combinatorial Efficient Global Optimization

###################################################################################
#' Build clustered model
#'
#' This function builds an ensemble of Gaussian Process model, where
#' each individual model is fitted to a partition of the parameter space.
#' Partitions are generated by.
#' 
#' @param x x
#' @param y y
#' @param distanceFunction distanceFunction
#' @param control control
#'
#' @return an object
#'
#' @export
###################################################################################
modelKrigingClust <- function(x, y, distanceFunction,control=list()){ 
  
  nsamples <- length(x) 
  con<-list(
    modelControl=list(), #controls passed to the core model
    minsize=30 #minimum number of observations for each partition to be modeled
  )
  con[names(control)] <- control
  control<-con
  
  #control$modelControl$target <- c("y","s") #required for weighting the ensemble predictions
  
  ## clustering by distance
  d <- distanceMatrix(x,distFun = distanceFunction)
  ## d_dist <- as.dist(distanceMatrix(x,distFun = distanceFunction))
  groups <- anticlust::balanced_clustering(d,max(floor(nsamples/control$minsize),2))
  unique_groups <- na.omit(unique(groups))
  
  ## can't get rid of the nasty NA, so just replace randomly
  nas <- is.na(groups)
  groups[nas] <- sample(unique_groups,sum(nas))

  ## build a model in each partition
  models <- list()
  for(i in 1:length(unique_groups)){
    selected <- groups == unique_groups[i]
    xi <- x[selected]
    yi <- y[selected]
    models[[i]] <-  modelKriging(x=xi,y=yi,distanceFunction=distanceFunction,
                                 control=control$modelControl
                                 )
  }
  
  ## ...
  fit <- list(fits=models)
  class(fit)<- "modelKrigingClust"
  return(fit)
}



###################################################################################
#' Clustered Kriging Prediction
#' 
#' Predict with a model fit resulting from \code{\link{modelKrigingClust}}.
#'
#' @param object fit of the clustered Kriging model ensemble (settings and parameters), of class \code{modelKrigingClust}.
#' @param newdata list of samples to be predicted
#' @param ... further arguments, currently not used
#'
#' @return list with function value (mean) \code{object$y} and uncertainty estimate \code{object$s} (standard deviation)\cr
#'
#' @seealso \code{\link{predict.modelKriging}}
#' @export
###################################################################################
predict.modelKrigingClust <- function(object,newdata,...){
  ## predict with each model
  predictions <- list()
  for(i in 1:length(object$fits)){
    fiti <- object$fits[[i]]
    fiti$predAll <- T
    predictions[[i]] <- predict(fiti,newdata)
  }
  
  ## convert predictions of s
  ps <- do.call(rbind,sapply(predictions,'[',"s"))
  
  ## convert predictions of y
  py <- do.call(rbind,sapply(predictions,'[',"y"))
  
  ## compute weights based on s
  psnegsquare <- ps^-2
  weights <- t(t(psnegsquare)/colSums(psnegsquare))
  
  ## compute ensemble prediction for s and y, weighted
  ensembley <- colSums(py*weights)
  ensembles <- sqrt(colSums(ps^2*weights^2))
  
  ## end
  list(y=ensembley,s=ensembles)
}