R/covariance.R

Defines functions updateCovModels getQLmodel fitSIRFk extractCovModels QLmodel fitCov doREMLfit fnGradREML fnREML reml setCovModel

Documented in fitCov fitSIRFk getQLmodel QLmodel reml setCovModel updateCovModels

# Copyright (C) 2017 Markus Baaske. All Rights Reserved.
# This code is published under the GPL (>=3).
#
# File: 	covariance.R
# Date:  	12/04/2017
# Author: 	Markus Baaske
#
# Define covariance structures: SIRF-k and Matern

#' @name
#'  setCovModel
#' 
#' @title 	
#' 	Set a covariance model
#' 
#' @description 
#' 	Set a covariance model for kriging the sample means of the involved statistics or for the variance matrix of the statistics.
#' 
#' @param model			name of covariance model
#' @param param			numeric vector, \code{NULL} (default), starting values of covariance parameters for estimation 
#' @param npoints		number of sample points already evaluated for covariance parameter estimation
#' @param var.sim		numeric vector, \code{NULL} (default), local simulation variances (so-called local nugget variances)
#' @param nugget		starting value for (global) nugget variance estimation
#' @param trend			trend order ID: either linear (=1) or quadratic (=2) for polynomial trend terms 					    
#' @param fixed.param	vector of names, corresponding to `\code{param}` of covariance parameters, which are hold fixed for covariance
#' 						parameter estimation
#' @param lower			lower bounds of covariance parameters for REML estimation 
#' @param upper			upper bounds of covariance parameters for REML estimation
#' @param ...			additional arguments which can be stored
#' 
#' @return Object of class \code{covModel}, a list of following elements
#'  \item{model}{ integer of covariance function}
#'  \item{param}{ estimtated covariance parameters}
#' 	\item{start}{ start point of REML estimation of covariance parameters}
#'  \item{trend}{ trend order}
#'  \item{fix.nugget}{ vector of (fixed) values used as local nugget variances}
#'  \item{free}{ index of free parameters for REML estimation}
#' 	\item{lower}{ lower bounds of covariance parameters for REML estimation}
#'  \item{upper}{ upper bounds of covariance parameters for REML estimation}
#'  \item{...}{ additional objects which can be stored}
#' 
#' @details The function defines a covariance model for kriging the sample mean values of an involved statistic. The covariance model
#'  (including a polynomial trend) defines the spatial dependence between different locations (points) of the
#'  parameter space. Currently, the function provides the generalized covariance models
#'  (`\code{sirfk}`, see \code{\link{fitSIRFk}}) of order \eqn{k=1,2} 
#'  and the Mat\eqn{\textrm{\'{e}}}rn covariance model with scale (or sill) parameter `\code{scale}`, smoothness parameter
#'  `\code{alpha}`, respectively, `\code{nu}`, and the range parameter `\code{rho}` defined only for the latter. A power exponential
#'  covariance model named "\code{powexp}" is also supported. 
#'  
#'  \subsection{Use of simulation variance}{
#'  If a vector of simulation variances is statically set by `\code{var.sim}` for each location these are used as (local)
#'  nugget variance estimations which account for the sampling variability due to the repeated measurements of the statistics
#'  by simulation. The length should match the number of locations `\code{npoints}` otherwise the given vector components
#'  are recycled. A global scalar valued nugget, which can captures the variance of the underlying random function, can be set by
#'  `\code{nugget}` as a starting value for the REML estimation procedure. Clearly, both types of nugget variances have a direct
#'  influence on the REML estimates.
#'  } 
#' 
#'  \subsection{Default parameters}{ 
#'  The default starting parameters are set to \code{("scale"=0.001,"alpha"=1.5)} for the `\code{sirfk}` model. The
#'  Mat\eqn{\textrm{\'{e}}}rn model uses the following parameters \code{("scale"=1.0,"nu"=2.5,"rho"=3.0)}. The default
#'  parameters for the power exponential covariance model are \code{"scale"=1.0}, (isotropic) range parameter
#'  \code{"phi"=1.0} and power \code{"kappa"=1.5} with \eqn{0<\kappa\leq 2}. 
#'  The corresponding lower and upper bounds are chosen such that the underlying random function remains
#'  twice continuously differentiable. Further, setting the names of the covariance parameters in `\code{fixed.param}`,
#'  excludes these parameters from subsequent REML estimations such that they remain unchanged. 
#' 
#'  The above settings are applicable for a wide range of statistics but, however, generally depend on the kind of statistics to be interpolated
#'  and thus have to be chosen carefully. Note that a valid (generalized) covariance model for kriging requires at least \eqn{q+2} design points
#'  for the trend order \eqn{k=1} and \eqn{1+(q+1)(q+2)/2} for \eqn{k=2} where \eqn{q} is the dimension of the unknown model parameter.
#'  }   
#'  
#' @examples 
#'  # set the standards sirf-2 covariance model
#'  setCovModel("sirfk",npoints=12)
#' 
#' @author M. Baaske
#' @rdname setCovModel
#' @export
setCovModel <- function(model = "sirfk", param = NULL, npoints = 0, var.sim = NULL,
						  nugget = 1.5e-4, trend = 2, fixed.param = NULL, 
						   lower = NULL, upper = NULL,...)
{	
	stopifnot(npoints>0)
	trend.nr <- pmatch(trend,c(1,2))
	if(anyNA(trend.nr)) {
	  msg <- paste0("Invalid trend order. Either choose `trend` as 1 (linear) or 2 (quadratic).")
	  message(msg)
	  stop(msg)
  	}
 	
    cov.list <- c("sirfk","matern","powexp","exp")
 	cov.nr <- pmatch(model, cov.list)
	if (is.na(cov.nr)) {
		msg <- paste("Unknown covariance model!",
				"Possible values are", paste(cov.list,collapse = ","))
		message(msg)
		stop(msg)
	}		
	
	fix.nugget <- NULL	
	if(!is.null(var.sim)) {		
	  if(!is.numeric(var.sim) || is.matrix(var.sim))
		 stop("`var.sim` has to be a numeric vector.")
	  fix.nugget <-
		# treated as local/global nugget
		if(length(var.sim) != npoints)
			as.numeric(unlist(lapply(list(var.sim), rep, length.out = npoints)))			
		else var.sim	 				
	}
	
	switch(model,
		"sirfk" = {			
			i <- min(which(trend < c(1,2,3), arr.ind=TRUE))
			param <-
			 if(is.null(param))
			  c("scale"=0.001,"alpha"=1.5,"nugget"=nugget)
			 else {			   
			   if(anyNA(match(names(param),c("scale","alpha"))))
				 stop("Invalid parameter vector for model `sirfk`.")
			   p <- c("scale"=0.001,"alpha"=1.5)  
			   p[names(param)] <- param 
			   c(p,"nugget"=nugget)
			 }
		 	
	        if(is.null(lower)) {
			  lower <- rep(-Inf,length(param))	
			}
			lower <-
			  if(anyNA(lower) | any(!is.finite(lower)))
				  c(1e-6,1.0,1e-4) else lower
			upper <-
			  if(is.null(upper)) {			  
				 c(5,i-1e-4,3)						
			  } else c(upper[1],min(upper[2],i-1e-4),upper[3])			 				
		},		
		"matern" = {
			param <-
			 if(is.null(param))
				c("scale"=1.0,"nu"=2.5,"rho"=3.0,"nugget"=nugget)
			 else {
				 if(anyNA(match(names(param),c("scale","nu","rho"))))
					 stop("Invalid parameter vector for model `matern`.")
				 p <- c("scale"=1.0,"nu"=2.5,"rho"=3.0)  
				 p[names(param)] <- param
			     c(p,"nugget"=nugget)
			 }
			if(is.null(lower) || is.null(upper)) {
				lower <- c(1e-6,2+1e-10,0.1,1e-4)
				upper <- c(10,3,10,10)
			}			
		},
		"powexp" = {
			param <-
			 if(is.null(param))
				param <- c("scale"=1.0,"phi"=1.0,"kappa"=1.5,"nugget"=nugget)
			 else {
				 if(anyNA(match(names(param),c("scale","phi","kappa"))))
					 stop("Invalid parameter vector for model `powexp`.")
				 p <- c("scale"=1.0,"phi"=1.0,"kappa"=1.5)  
				 p[names(param)] <- param				 
				 c(p,"nugget"=nugget)
			 }
	 
			if(is.null(lower) || is.null(upper)) {
				lower <- c(1e-6,1e-6,1e-4,1e-4)
				upper <- c(3,3,1.9999,10)
			}			
		},
		"exp" = {
			param <-
			if(is.null(param))
				param <- c("scale"=1.0,"phi"=1.0,"nugget"=nugget)
			else {
				if(anyNA(match(names(param),c("scale","phi"))))
					stop("Invalid parameter vector for model `exp`.")
				p <- c("scale"=1.0,"phi"=1.0)  
				p[names(param)] <- param				 
				c(p,"nugget"=nugget)
			}
			
			if(is.null(lower) || is.null(upper)) {
				lower <- c(1e-6,1e-6,1e-4)
				upper <- c(3,3,10)
			}			
		}
	)	
	if(length(param) != length(lower) ||
	   length(param) != length(upper)) {
	  stop("`lower` or `upper` bound lengths must be equal to the length of the covariance parameter vector.")
	}
	start <- param
	if(any(start<lower) || any(start>upper))
		stop("At least one parameter does not match constraints. Check `lower` and `upper`.")
	
			
	free <-  
	 if(is.null(fixed.param)){
		 seq(length(start))
	 } else {
		id <- which(is.na(match(names(start),fixed.param)))
		if(length(id) == 0) {
		  cat("All parameters are excluded estimation.")
		} else {
		  start <- start[id]	    # new start
		  lower <- lower[id]		# new bounds
		  upper <- upper[id]		  
		} 
		id
	 }  	
		
	structure(list("model"= cov.nr,
				   "param"= param,
				   "start"= start,				  				   				   
				   "trend"= trend,
				   "fix.nugget"= fix.nugget,	
				   "free"= free,
				   "lower"= lower,
				   "upper"= upper,...),
			class="covModel"
	)
	
}

#' @name 		reml
#' 
#' @title 		Restricted maximum likelihood (REML)
#' 
#' @description Calculate the REML value (witout constant term) for a given covariance model and data 
#' 
#' @param models   	 object of class \code{krige} (list of covariance models) or class
#'					 	 \code{covModel} (a single covariance model), see \code{\link{setCovModel}}
#' @param pars 	     covariance parameter vector (including global scalar nugget value)
#' @param data	  	 data frame of simulated statistics, each column corresponds to a 
#' 					 single covariance model in `\code{models}`
#' @param Xs	  	 matrix of sample points
#' @param verbose	 logical, if \code{TRUE}, print intermediate output
#' 
#' @return List of REML function values.
#' 
#' @details Given a list of covariance models the function calculates the REML function values at
#'   the covariance parameter `\code{pars}`.
#' 
#' @examples
#' 	data(normal)
#' 
#'  # extract the sample points (the design)
#'  X <- as.matrix(qsd$qldata[1:2])
#' 
#'  # get the observed statistic
#'  T <- qsd$qldata[c("mean.T1")]
#'  reml(qsd$covT[1],pars=c(1e-4,1.5,0.1),T,X)
#' 
#' @author M. Baaske
#' @rdname reml
#' @export
reml <- function(models, pars, data, Xs, verbose = FALSE) {	
	stopifnot(is.data.frame(data))
	stopifnot(length(models) == ncol(data))
	
	lapply(1:length(models),
		function(i) {
			tryCatch({		
				 F <- .Call(C_Fmat,Xs,models[[i]]$trend)	
				 P <- .Call(C_Pmat,F)	
				 y <- crossprod(P,data[[i]]) # t(P)%*%z
				 fnREML(pars,y,Xs,P,models[[i]],verbose=verbose)
				}
			  ,error = function(e) {
				  msg <- .makeMessage("REML function evaluation failed: ",conditionMessage(e))
			  	  message(msg) 
				  structure( list( message=msg,  call = sys.call() ), class=c("error","condition"), error=e )
			 	}	  
			)			
		})
}

# intern, covModel has to be initialized! 
fnREML <- function(p, y, Xs, P, model, free = seq(length(p)), verbose = FALSE)
{	
	 tryCatch({			
		model$param[free] <- p		
		# covariance matrix
		Cmat <- .Call(C_covMatrix,Xs,model)		
		if (!is.matrix(Cmat) || anyNA(Cmat) || !is.numeric(Cmat))
		  stop("Error in covariance matrix: non-numeric values.")								
			
	 	msg <- NULL	
	 	W <- crossprod(P, Cmat %*% P)
		rc <- rcond(W)
				
	    if(rc > 1e-3){
			Wc <- try(chol(W),silent=TRUE)
			if(!inherits(Wc,"try-error") && all(diag(Wc)>0) ){
				w  <- backsolve(Wc,y,transpose=TRUE)	
				if(!inherits(w,"try-error")){
				  return (
				   structure(
					as.numeric(.5*(sum(w^2)) + sum(log(diag(Wc)))),  # beware of brackets: 0.5*2*sum(log(diag(Wc)))
				  	 "info"=list("rcond"=rc,"msg"=msg, "p"=p)) 
	              )
			  	} else {
					if(verbose) 
					 warning("Could not backsolve in `fnREML`.")
				}
			} 
		} else {		  
			msg <- paste0(c("reciprocal condition number (", rc,") of projected covariance matrix `W` is near zero at covariance parameter \n\t ",
					format(p, digits=6, justify="right"),"\n"," which might be unreliable."), collapse = " ")
			if(verbose)
				warning(msg)
		}				
		z <- try(gsiSolve(W,y,use.solve=FALSE),silent=TRUE)
		if(inherits(z,"try-error"))
		  stop("`gsiSolve` failed. Cannot continue REML estimation.")
	  	detW <- try(det(W),silent=TRUE)
		if(inherits(detW,"try-error") || detW < 0)
		 stop(.makeMessage("Could not compute determinant: ",detW," or negative value. \n"))
	 	if(verbose && detW < 1e-17)
		 warning("Determinant of projection matrix `W` (REML) is near zero.\n")
		
		return(	
		    structure(
				as.numeric( .5*(t(y) %*% z + log(detW)) ),
				 "info"=list("rcond"=rc,"msg"=msg, "p"=p) )	
		)		
	  
	} ,error = function(e) {
		 stop(.makeMessage("Error in function 'fnREML': ",	conditionMessage(e)))		  	
		}
	)		
	
}  

## TODO: Numerical derivatives of log likelihood
##  -> use gradient of loglik
## REML ->alpha
## alpha <- 1.5
## phi <- 1.0
## D <- as.matrix(dist(X))
## diag(D)<- (fix.nugget+nugget)
## H <- log(phi*D)
## exp(2*alpha*H)*2*H
#' @importFrom nloptr nl.grad
fnGradREML <- function(p, y, Xs, P, model, free = NULL, verbose = FALSE) {
	list("objective"=fnREML(p,y,Xs,P,model,free,verbose),
		 "gradient"=nloptr::nl.grad(p, fnREML, heps = .Machine$double.eps^(1/3),y,Xs,P,model,free)) 
}  

## TODO add data as parameter
# arguments '...' manipulate global options for nloptr 
#' @importFrom nloptr nl.opts
doREMLfit <- function(model, Xs, opts, verbose = FALSE )
{
	# return if all parameters are fixed
	if(!is.null(model$free) && length(model$free)==0L) {
	  return(
		structure(
		  list(model = model,convergence = 1L),
		optres=NULL, class = "reml") )
	}	
	err <- NULL	
	fn <- fnREML
	if(!is.null(opts$algorithm) &&
	   opts$algorithm == "NLOPT_LD_LBFGS") {
	   message("Caution: using `LBFGS` in REML function might not be stable.")
	   fn <- fnGradREML
 	}
	 
 	tryCatch({		
        nms <- names(model$param)				
		Fmat <- .Call(C_Fmat,Xs,model$trend)		
		P <- .Call(C_Pmat,Fmat)	
		
		# transform to mean zero data				
		y <- crossprod(P,as.numeric(model$dataT))
		
		model$dataT <- NULL		
		p0 <- .PROJMED(model$start,model$lower,model$upper)
		
		res <- nloptr::nloptr(p0, fn, lb = model$lower, ub = model$upper, opts = opts,
						y = y, Xs = Xs, P = P, model = model, free = model$free,
						 verbose = verbose)
		msg <- "Normal convergence."
		if(inherits(res,"error") || is.null(res) || anyNA(res$solution)){
			msg <- .makeMessage("Function call to 'nloptr' failed.")				
			message(msg)
			return(.qleError(message=msg,call=match.call(),error=res))
		}
		# do a final local search
		if(!is.null(opts$local_opts)){
			if(length(opts$local_opts) > 0L) {		
				locopts <- nloptr::nl.opts()
				locopts[names(opts$local_opts)] <- opts$local_opts 
			} else {
				locopts <- list("algorithm" = "NLOPT_LN_COBYLA","ftol_rel" = 1.0e-7,
								"xtol_rel" = 1.0e-6, "maxeval" = 100)
			}	
			if(verbose)
			  message("Do a final local search of covariance parameters.")
			res0 <- nloptr::nloptr(res$solution, fn, lb = model$lower, ub = model$upper, opts = locopts,
					 y = y, Xs = Xs, P = P, model = model, free = model$free,
					  verbose = verbose)
			if(inherits(res0,"error") || is.null(res0) || anyNA(res0$solution)){
			   warning(.makeMessage("Local function call to 'nloptr' failed after global optimization."))
			   res$final <- res0
		    } else res <- res0
		}
		
		converged <- FALSE
		sol <- res$solution			
		if(res$status >= 0L) {
			converged <- TRUE
			model$param[model$free] <- sol						   			   
	    } else {
		   verbose <- TRUE
		   msg <- .makeMessage("Estimation of covariance parameters did not converge.")		   
		   message(msg)
	    }		
	    structure(
		    list(model=model,
				 convergence=converged,
				 message=msg),
		  optres = if(verbose) res else NULL, class = "reml") 
				 
	 }, error = function(e){
			 msg <- .makeMessage("Nloptr error fitting covariance parameters: ",
					  conditionMessage(e))
			 message(msg)
			 structure(
				list(model=model,
					 convergence=FALSE,
					 message=msg,
					 call=sys.call()),
			  error=e)			  		
		}
	) 	 
}


#' @name fitCov
#' 
#' @title Fitting covariance models by REML estimation
#' 
#' @description The function estimates the (hyper)parameters of the covariance models by
#' 	  the \emph{Restricted Maximum Likelihood} (REML) method.
#' 
#' @param models  	 object either of class \code{krige}, a list of covariance models, or of
#' 	 				 class \code{covModel}, that is, a single covariance model
#' @param Xs	 	 matrix of sample points, the design
#' @param data		 data frame of simulated sample means of statistics
#' 					 first column corrsponds to the first model in the list `\code{models}` and so forth
#' @param controls	 list of control parameters, see \code{\link[nloptr]{nloptr}}
#' @param cl		 cluster object, \code{NULL} (default), of class "\code{MPIcluster}", "\code{SOCKcluster}", "\code{cluster}"
#' @param verbose 	 logical, \code{TRUE} for intermediate output
#' 
#' @return An object of class \code{reml} which consists of a list of named lists
#'  (`\code{model}`, `\code{convergence}`) each storing a fitted covariance model itself
#'  together with the optimization results from \code{\link[nloptr]{nloptr}} as an attribute
#'  named `\code{optres}`. The default method for estimating the covariance parameters is
#'  \code{\link[nloptr]{mlsl}} which uses random starting points and thus could produce different results if
#'  it is run multiple times. If the results strongly vary, then the corresponding REML function might have many
#'  local minima which precludes the use of this default algorithm and another one, e.g. `\code{NLOPT_GN_DIRECT}`
#'  (see \code{\link[nloptr]{nloptr.print.options}}), might lead to better results. 
#' 
#' @details The function fits a list of covariance models using the REML method. In order to avoid singularities
#'  of the so-called trend matrices make sure to use at least the minimum required number of sample points stored in
#'  `\code{Xs}` which depends on the defined trend order, see \code{\link{setCovModel}}.
#' 
#' @examples 
#' data(normal)  
#' 
#' # fit 1st statistic and get REML results
#' fitCov(qsd$covT[1],
#'        Xs=as.matrix(qsd$qldata[1:2]),
#'        data=qsd$qldata["mean.T1"],verbose=TRUE)
#'   
#' @seealso \code{\link{setCovModel}} 
#' 
#' @author M. Baaske
#' @rdname fitCov 
#' @export 
fitCov <- function(models, Xs, data, controls = list(),
			     	  cl = NULL, verbose = FALSE) {
		
	if(!is.data.frame(data))
		stop("Expected argument `data` of class `data.frame`.")	
	if(!is.matrix(Xs))
		stop("Expected argument `Xs` to be  a matrix of sample locations.")
			
	if(length(controls)>0L) {		
		opts <- nloptr::nl.opts()
		opts[names(controls)] <- controls
	} else {
		opts <- list("algorithm" = "NLOPT_GN_MLSL",
				"local_opts" = list("algorithm" = "NLOPT_LN_COBYLA","ftol_rel" = 1.0e-6,
						"xtol_rel" = 1.0e-6,"maxeval" = 1000),
				"maxeval" = 200, "xtol_rel" = 1.0e-6, "ftol_rel" = 1.0e-6, "population"=0)	
	}	
	for(i in 1:length(models))
	 models[[i]]$dataT <- as.numeric(data[[i]])
		
 	mods <- doInParallel(models, doREMLfit, Xs=Xs, opts = opts,
			  cl=cl, verbose=verbose)
	  
	if(inherits(mods,"error")) {
		msg <- paste0("REML estimation failed: ",conditionMessage(mods),"\n")
		message(msg)
		return(.qleError(message=msg,
				call=match.call(),error=mods))
	}	
	errId <- which(sapply(mods,function(x) .isError(x)))
	if(verbose) {	  
	  if(any(errId))
		message(paste(c("Failed fitting covariance models with index: ",as.character(errId)), collapse=" ")) 
	  else {
		id <- which(sapply(mods,function(x) x$convergence))
		if(!all(id)) {
		 message(paste(c("REML failed to converge: ",as.character(id)), collapse=" ")) 
		} else
		 message("Successfully fitted covariance parameters.","\n")		
	  }
	}	
	structure(mods,
		opts = opts,
		error = if(length(errId) > 0L) errId else NULL,
		class = "QLFit")
}

#' @name QLmodel
#' 
#' @title Construct quasi-likelihood approximation 
#' 
#' @description Aggregate and construct the data for quasi-likelihood estimation
#' 
#' @param qldata		data frame of (initial) simulation results (see \code{\link{setQLdata}})
#' @param lb		    numeric vector of lower bounds defining the (hyper)box
#' @param ub 			numeric vector of upper bounds defining the (hyper)box
#' @param obs	    	numeric vector of observed statistics
#' @param mods			list of (fitted) covariance models (see \code{\link{fitSIRFk}}) 
#' @param nfit			number of cycles, \code{nfit=1} (default), after which covariance
#' 						parameters are re-estimated otherwise re-used 
#' @param cv.fit 		logical, \code{TRUE} (default), whether to re-fit CV models (re-estimate covariance parameters)	
#' @param var.type  	name of the variance approximation method (see \code{\link{covarTx}})
#' @param useVar    	logical, \code{TRUE} (default), whether to use prediction variances (see details)
#' @param criterion 	global criterion function for sampling and minimization, either "\code{qle}" or "\code{mahal}"				    	
#' @param verbose       logical, \code{FALSE} (default), whether to give further output 
#' 
#' @return An object of class \code{\link{QLmodel}} which stores the data frame of simulation results, bounds on
#'  the parameter space, covariance models for kriging, vector of observed statistics as well as options for
#'  kriging and fitting. 
#' 
#' @details The function aggregates all required information for quasi-likelihood estimation, stores the fitted
#'   covariance models of the sample means of the statistics and the type of variance matrix approximation. For an advanced setup
#'   of the estimation procedure and more involved statistical models this function explicitly offers the data structure to construct
#'   individual covariance models for each statistic as defined by \code{\link{setCovModel}}. The user has the choice whether or not to
#'   make use of of kriging prediction variances by `\code{useVar}` to account for the simulation error when constructing
#'   the approximation of the variance matrix and the quasi-score function. If \code{TRUE}, then a kriging procedure calculating
#'   prediction variances is automatically used. Otherwise the so-called \emph{dual} approach is employed which has some computational
#'   advantage if prediction variances are not required. 
#' 
#' @examples 
#' 
#' data(normal)
#' 
#' # We simply re-use the stored normal data and fit again: 
#' # fit generalized covariance model to the data using
#' # simulation variances as local nugget variances
#' mods <- fitSIRFk(qsd$qldata, verbose=TRUE)
#' 
#' # construct QL approximation model
#' qsd <- QLmodel(qsd$qldata,qsd$lower,qsd$upper,
#' 			    c("T1"=2,"T2"=1),mods)
#' 
#' @seealso \code{\link{getQLmodel}}, \code{\link{updateCovModels}} 
#' 
#' @author M. Baaske
#' @rdname QLmodel
#' @export 
QLmodel <- function(qldata, lb, ub, obs, mods, nfit = 1, cv.fit = TRUE,
		    var.type = c("wcholMean","cholMean","wlogMean","logMean","kriging","const"),
				useVar = TRUE, criterion = c("qle","mahal"), verbose = FALSE)
{	
	if(!inherits(qldata,"QLdata"))
	 stop("expected argument `qldata` of class `QLdata`.")
	if(missing(lb) || missing(ub))
	 stop("Arguments `lb` and `ub` are missing.")
 	
 	dx <- attr(qldata,"xdim")
    if(!is.numeric(lb) || !is.numeric(ub) ||
	   length(lb)!=length(ub) || dx != length(ub))
  	  stop("Dimensions of `lb` or `ub` do not match.")
 	
 	obs <- unlist(obs)
	if(anyNA(obs) | any(!is.finite(obs)))
	 stop("`NA`,`NaN` or `Inf`values detected in argument `obs.")
	if(!is.numeric(obs))
	  stop("Argument `obs` must be a (named) numeric vector or list.")
  	stopifnot(!is.null(mods$covT))
	stopifnot(class(mods)=="QLFit")
	
	if(length(mods$covT) != length(obs))
	  stop("Number of covariance models `covT` and length of observations vector `obs` must equal.")
  	var.type <- match.arg(var.type)
	criterion <- match.arg(criterion)
		
	if(is.null(mods$covL) && var.type == "kriging")
	  stop("Covariance models for variance matrix interpolation must be set for argument \'var.type\'.")
	if(!is.numeric(nfit) || length(nfit)>1L )
	 stop("Argument 'nfit must be numeric of length one.")	
	
 	covT <- .extractCovModels(mods$covT,verbose)
	stopifnot(class(covT)=="krige")
	
	covL <- NULL
	if(!is.null(mods$covL)){		
		covL <- .extractCovModels(mods$covL,verbose)
		stopifnot(class(covL)=="krige")
	}
	# reml optimization options 
	opts <- attr(mods,"opts")
	if(is.null(opts) || length(opts) == 0L){
		opts <- list("algorithm" = "NLOPT_GN_MLSL",
				"local_opts" = list("algorithm" = "NLOPT_LN_COBYLA","ftol_rel" = 1.0e-6,
						"xtol_rel" = 1.0e-6,"maxeval" = 1000),
				"maxeval" = 200, "xtol_rel" = 1.0e-6, "ftol_rel" = 1.0e-6, "population"=0)			  
	}
	# minimum required sample size
	minN <- ifelse(min(sapply(covT,	function(x) x$trend)) < 2, dx+2, (dx+1)*(dx+2)/2+1)
	if(nrow(qldata)<minN) {
	 stop(paste0("Choose the size of the initial sample for parameter dimension ",dx,
			" at least of size: ",minN))
 	}	
	
	structure(
	    list("qldata" = qldata,
			 "lower" = lb,
			 "upper" = ub,
			 "covT" = covT,
			 "covL" = covL,			
			 "obs" = obs,		 
			 "var.type" = var.type,			 
			 "krig.type" = if(useVar) "var" else "dual",
			 "criterion" = criterion,			 
			 "minN" = minN,
			 "nfit" = nfit,
			 "cv.fit"=cv.fit),
	  opts = opts,
	  class="QLmodel"
	)
}

# intern
.extractCovModels <- function(covs, verbose = FALSE) {
	if(is.null(covs))
	  return (NULL)
    # which one procudces errors in fitting?
	errId <- which(sapply(covs,function(x) .isError(x)))
	
	if(length(errId)>0L)
	  message(.makeMessage("A total of ",length(errId)," errors detected in fitted covariance models.")) 
	
    structure(
		lapply(covs,
		 function(x) {
			if(verbose)
			 structure(x$model,"optres"=attr(x,"optres"))
			else x$model
		 }
	    ), error = if(length(errId)>0L) errId else NULL,
	  class = "krige"
 	)	
}

#' @name fitSIRFk
#' 
#' @title Estimation of covariance parameters
#' 
#' @description Fit a generalized covariance model to simulation data
#' 
#' @param qldata		object of class \code{QLdata}, a data frame from \code{\link{setQLdata}}
#' @param set.var 		logical vector of length one or equal to the number of covariance models;
#' 						for values \code{TRUE} (default), set simulation variances as local nugget variances
#' 						for the corresponding covariance model/statistic 
#' @param var.type      name of variance matrix approximation type (see \code{\link{covarTx}})  
#' @param var.opts	    list of arguments passed to \code{\link{setCovModel}}
#' 						(only if `\code{var.type}`="\code{kriging}" and ignored otherwise)
#' @param intrinsic 	logical vector, \code{FALSE} (default), of length one or equal to the number of Cholesky
#' 					    decompositions of variance matrices; as default use an internal nugget variance estimate (see details)
#' 						for kriging the variance matrix of the involved statistics
#' @param ...			arguments passed to \code{\link{setCovModel}}
#' @param cl			cluster object, \code{NULL} (default), of class "\code{MPIcluster}", "\code{SOCKcluster}", "\code{cluster}"
#' @param controls		list of control parameters passed to \code{\link[nloptr]{nloptr}} for local minimization
#' @param verbose		if \code{TRUE}, print intermediate results
#' 
#' @return A list of fitted covariance models for kriging the sample means of statistics named `\code{covT}` and optionally
#'  the variance matrix of statistics, `\code{covL}`. The object also stores the reml optimization parameters `\code{controls}`. 
#' 
#' @details The function estimates the parameters of a covariance model using the REML method for kriging
#'   the sample means of the statistics and kriging the variance matrix of statistics unless `\code{var.type}`
#'   equals "\code{const}". By default it uses the covariance model derived from a (self-similar) intrinsic random function, that is,
#'   the `\code{sirfk}` function of order \eqn{k} (see, e.g. [1]) with \eqn{k=1,2}, for all statistics (including a default quadratic drift term
#'   \eqn{k=2}). The user can also define different covariance models for each statistic separately (see below). Other covariance models can be used by their
#' 	 name `\code{model}` which is passed to the function \code{\link{setCovModel}}. Kriging the variance matrix always uses the `\code{sirfk}` covariance model. 
#'    		
#'   Argument `\code{var.opts}` only sets the options for the covariance models for kriging the variance matrix if this is the users prefered
#'   type of approximation. Other optional arguments, e.g., `\code{var.sim}` for the statistics, `\code{var.opts$var.sim}` for kriging the variance matrix,
#'   specify the local or global  \dfn{nugget} values for each sample point depending on whether or not `\code{set.var}` (used for kriging the statistics)
#'   equals \code{TRUE}. Both are passed to \code{\link{setCovModel}} and must be data frames of lengths (number of columns) corresponding to the number of covariance
#'   models of statistics and, respectively, to the number of \emph{Cholesky} decomposed terms in case of kriging the variance matrix.
#'   If `\code{set.var}` equals \code{TRUE} (default), then local nugget variances are estimated by the variance of the sample average of the statistics.
#'   Otherwise the values given in `\code{var.sim}` are used as fixed `nugget` variances and replicated to match the number of sample points if required.
#'  
#'   The same principle applies in case of kriging the variance matrix. If `\code{intrinsic}` equals \code{TRUE}, then local nugget variances
#'   for each of the variance-covariances of the  of the statistics are estimated by a bootstrapping procedure. Otherwise the values given by `\code{var.opts$var.sim}`
#'   (of length one or equal to the number of corresponding sample points) are used directly as local estimates (which then must correspond to
#'   the other Cholesky decomposed terms). A global nugget value can be also estimated during the REML estimation which is the default option for both
#'   cases unless this parameter is excluded from the covariance parameter estimation (see \code{\link{setCovModel}}). The default optimization algorithm for
#'   estimating the covariance parameters is the random starting point method \code{\link[nloptr]{mlsl}} followed by a final local search by the same local algorithm.
#'   Note that in this case the estimated parameters may vary when starting the REML procedure several times since starting points are chosen as random. All
#'   options for the optimization can be modified by the argument `\code{controls}`.
#' 
#'   Note that the returned object can also be constructed manually and passed as an input argument to
#'   \code{\link{QLmodel}} in case the user prefers to set up each covariance model separately. In this case, first use \code{\link{setCovModel}} to construct
#' 	 the covariance model, then estimate the parameters by \code{\link{fitCov}} and pass a list of fitted covariance models to function \code{\link{QLmodel}}. 	 
#'  
#' 	 Please see function \code{\link{QLmodel}} for an example. 
#' 
#' @seealso \code{\link{setCovModel}}, \code{\link{fitCov}},  \code{\link{QLmodel}}
#' 
#' @author M. Baaske
#' @rdname fitSIRFk
#' @export
fitSIRFk <- function(qldata, set.var = TRUE, var.type = "wcholMean",
						var.opts = list("var.sim"=1e-6), intrinsic = FALSE, ...,
						 controls = list(), cl = NULL, verbose = FALSE)
{	
	args <- list(...)
	stopifnot(is.data.frame(qldata))
	stopifnot(is.logical(set.var) && is.logical(intrinsic))
	
	if(length(args)>0L) {		 
		nms <- formals(setCovModel)
		.checkfun(setCovModel, c(nms,args))			
 	}
	
	xdim <- attr(qldata,"xdim") 									# dimension of parameter to estimate!
	nsim <- attr(qldata,"nsim")										# number of simulation replications
	Xs <- data.matrix(qldata[seq(xdim)])	
	dataT <- qldata[grep("^mean.",names(qldata))]					# simulated statistic data
		
	np <- nrow(Xs)	
	nstat <- ncol(dataT)	
	
	# all cov models are equal	
	useVarSim <- !is.null(args$var.sim)	
	# if not used as a fixed nugget: set.var == FALSE	
	set.var <- rep(set.var,length.out=ncol(dataT))
	dfvar <-
	 if(useVarSim) {
	  if(anyNA(args$var.sim))
		 stop("local nugget variance vector has `Na`s for kriging statistics.")
	  rep(as.data.frame(args$var.sim),length.out=ncol(dataT))		
	} else NULL
	
	covT <- 
		lapply(1:ncol(dataT),
		     function(i){
			   args$var.sim <-
				 if(set.var[i]) {					  
					  qldata[[xdim+nstat+i]]/nsim
				 } else if(useVarSim && any(dfvar[[i]]>0)) {				 
					  dfvar[[i]]										# numeric vector of length equal number of locations
				  } else NULL
				 fnargs <- c(list("dataT"=dataT[[i]],	      		    # temporarly add the data				  				  			  
								  "npoints"=np,
								  "type"="statcov"), args)
				 do.call(setCovModel,fnargs)					   
			}
		)

	 covL <- NULL	 
	 if(var.type == "kriging"){
		 # check input
		 args <- var.opts
		 if(length(args)>0L) {
			 nms <- formals(setCovModel)			
			.checkfun(setCovModel, c(nms,args))				 
	 	 }		 	   		   		 
  		 useVarSim <- !is.null(args$var.sim)
		 Lvec <- qldata[grep("^L+",names(qldata))]		 
		 
		 # individually set intrinsic noise terms as local nugget variances
		 # for each covariance model of Cholesky decomposed terms
		 intrinsic <- rep(intrinsic,length.out = ncol(Lvec))		 
		 dfvar <- 
		   if(useVarSim) {
			   if(anyNA(args$var.sim))
				   stop("local nugget variance vector has `Na`s for kriging variance matrix.")
			   rep(as.data.frame(args$var.sim),length.out=ncol(Lvec))			   
		   } else { NULL }
		 
   		 # find number of additional columns in `qldata`
   		 M <- if(attr(qldata,"Nb")>0) (nstat*(nstat+1))/2 else 0		 
		 # first M columns are Cholesky terms,
	     # 2nd are bootstrap variances
		 covL <- lapply(1:(ncol(Lvec)-M),
				  function(i)  {
				     args$var.sim <-
					 if(intrinsic[i] && M>0) {						 							 
						 if(any(Lvec[[i+M]] < 0) || anyNA(Lvec[[i+M]])){
							if(verbose)
							  message("Bootstrap variance is negative or has `Na`s. So we set a default nugget variance value.")
						    # set small value anyway
							as.numeric(dfvar[[i]]) 
						 } else {
							 # square root of nugget variance as this corresponds
							 # to a single value of the Cholseky decomposition 
							 sqrt(Lvec[[i+M]])
						 }
					 } else if(useVarSim && any(dfvar[[i]]>0)) {						 
						 as.numeric(dfvar[[i]]) 
					 } else NULL
					 	 	    
					 fnargs <- c(list("dataT"=Lvec[[i]],									  		  
									  "npoints"=np,
									  "type"="kriging"), args)							   	  
					 do.call(setCovModel,fnargs)					 
				  }
		 )	 
	 }		 	 
	 # (default) reml optimization options 
	 if(length(controls) > 0L) {		
		 opts <- nloptr::nl.opts()
		 opts[names(controls)] <- controls
	 } else {
		 opts <- list("algorithm" = "NLOPT_GN_MLSL",
				  "local_opts" = list("algorithm" = "NLOPT_LN_COBYLA","ftol_rel" = 1.0e-6,
						 "xtol_rel" = 1.0e-6,"maxeval" = 1000),
				  "maxeval" = 200, "xtol_rel" = 1.0e-6, "ftol_rel" = 1.0e-6, "population"=0)		  
	 }
	 	 
	 # REML fit covariance models (statistics and variance matrices)
	 mods <- doInParallel(c(covT,covL), doREMLfit, Xs=Xs, opts=opts,
			 	cl=cl, verbose=verbose)
		
	 if(inherits(mods,"error")) {
		msg <- paste0("REML estimation failed: ",conditionMessage(mods),"\n")		
		message(msg)
		return(.qleError(message=msg,
				 call=match.call(),error=mods))
	 }
	 errId <- which(sapply(mods,function(x) .isError(x)))
	 if(any(errId)) {
		 msg <- paste(c("Failed fitting covariance parameters: ",
				   as.character(errId)), collapse=" ")
   		 message(msg)
		 return(.qleError(message=msg,error=mods[errId]))
	 } else {
		 if(verbose)
		   cat("Successfully fitted covariance parameters.\n")		 
	 }
	ret <- structure(
			  list("covT" = mods[1:nstat],
			  	   "var.type" = var.type),
		     opts = opts,
			 error = if(length(errId)>0L) errId else NULL,
			 class = "QLFit")	
	 
	if(!is.null(covL))
	  ret$covL <- mods[(nstat+1):length(mods)]
	  	
	return ( ret )	
}

#' @name getQLmodel
#' 
#' @title Setup the quasi-likelihood estimation model  
#' 
#' @description  Setup the quasi-likelihood model data
#'
#' @param runs   	object of class \code{simQL} as simulation results from \code{\link{simQLdata}}
#' @param lb		lower bounds defining the (hyper)box
#' @param ub 		upper bounds defining the (hyper)box
#' @param obs		numeric vector of observed statistics
#' @param X   		matrix of sample locations (model parameters)
#' @param useVar   	logical, \code{TRUE} (default), whether to use prediction variances
#' @param criterion the criterion function to be minimized for parameter estimation (see \code{\link{qle}})
#' @param ...		arguments passed to \code{\link{fitSIRFk}}, \code{\link{setQLdata}}, \code{\link{setCovModel}} 
#'  					and \code{\link{QLmodel}} for fitting kriging covariance models
#'  
#' @return Object of class \code{\link{QLmodel}}
#' 
#' @details The function is a wrapper to \code{\link{setQLdata}} and \code{\link{fitSIRFk}}
#'  in order to setup the data required for estimating the model parameters.
#' 
#' @examples
#' 
#' data(normal)
#' 
#' # simulate model at a minimum of required design points
#' sim <- simQLdata(sim=qsd$simfn,nsim=5,N=8,
#' 			 method="maximinLHS",lb=qsd$lower,ub=qsd$upper)
#' 	 
#' # true and error-free observation
#' obs <- structure(c("T1"=2,"T2"=1), class="simQL")
#' 
#' # construct QL approximation model
#' qsd <- getQLmodel(sim,qsd$lower,qsd$upper,obs,var.type="wcholMean")
#' 
#' 
#' @seealso \code{\link{simQLdata}}, \code{\link{QLmodel}}, \code{\link{fitSIRFk}} 
#' 
#' @author M. Baaske
#' @rdname getQLmodel
#' @export
getQLmodel <- function(runs, lb, ub, obs, X = NULL, useVar = TRUE, criterion = "qle", ...)
{	
	args <- list(...)
	verbose <- isTRUE(args$verbose)
		
	tryCatch({
        if(.isError(runs))
		  stop("Simulations have errors. Please check the input argument `runs`.")
		if(verbose)
		  cat("Collect data for fitting covariance models of statistics.\n")
	  		  
	  	id <- which(is.na(pmatch(names(args),names(formals(setQLdata)))))
		args.tmp <-
		 if(length(id)>0L){
		   if( (c("Nb") %in% names(args)) && !isTRUE(args$intrinsic) )	 
		     args$Nb <- 0  # no bootstrap anyway if not "intrinsic" equals TRUE
		   args[-id]
	     } else NULL 
 		# construct all data
	    qldata <- do.call(setQLdata,c(list(runs,X),args.tmp))
		if(.isError(qldata))
			return(qldata)
		
		# fitting statistics
		if(verbose)
		 cat("Fitting covariance models...\n")	 	
	    id <- which(is.na(pmatch(names(args),names(c(formals(fitSIRFk),formals(setCovModel))))))
	 	args.tmp <-
		 if(length(id)>0L) {
			args[-id]
		 } else args	    
		
		# fitting			 	
	    mods <- do.call(fitSIRFk,c(list(qldata),args.tmp))	
		if(.isError(mods)) {
			message(.makeMessage("Failed fitting covariance models: ","\n"))			
			return(mods)
		}		  
		if(verbose)
		  cat("Setup QL approximation model...\n")
	  	id <- which(is.na(pmatch(names(args),names(formals(QLmodel)))))
	  	if(length(id) > 0L) {
		  args <- args[-id]		  
	 	}
	    do.call(QLmodel,
			 c(list(qldata,
					lb,ub,
			  	    obs,
				    mods,				    
					useVar=useVar,
				    criterion=criterion),
		     args))
		
	}, error = function(e) {
		msg <- .makeMessage("Failed to setup QL model: ",conditionMessage(e))
		message(msg)
		return(.qleError(message=msg,
				 call=match.call(),error=e))	
	   }
	)
}

#' @name updateCovModels
#' 
#' @title Update covariance models
#' 
#' @description The function updates the current covariance models
#'  stored in `\code{qsd}`.
#' 
#' @param qsd			object of class \code{\link{QLmodel}} which is to be updated
#' @param nextData		object of class \code{QLdata} which includes new simulation results
#' @param fit 			logical, if \code{TRUE} (default), re-estimate covariance parameters
#' @param cl			cluster object, \code{NULL} (default), of class "\code{MPIcluster}", "\code{SOCKcluster}", "\code{cluster}"
#' @param controls	    list of control parameters passed to \code{\link[nloptr]{nloptr}}
#' @param verbose 		logical, whether to print intermediate information
#' 
#' @return Object of class \code{\link{QLmodel}} as a list of updated covariance models
#' 
#' @details The function updates both, the covariance models for kriging the statistics, and, if applicable,
#'  the ones for kriging the variance matrix of statistics based on the new data \code{nextData}. In practice, the user hardly
#'  needs to call this function except for empirical studies of how additional sample points might influence the overall predictive
#'  quality of the quasi-score and/or criterion function approximations.
#' 
#'  If `\code{fit}` equals \code{TRUE}, then the function re-estimates the covariance parameters for each statistic separately
#'  each time a total of `\code{qsd$nfit}` new sample points have been added. Thus, we can choose whether to fit the updated
#'  covariance models (by the REML estimation) each time, e.g. during the estimation by \code{\link{qle}} if `\code{qsd$nfit}`=1, or after
#'  each 2nd, 3rd, and so on newly added point in order to limit the computational overhead.
#' 
#' @examples
#' 
#' data(normal)
#' 
#' # old design
#' X <- as.matrix(qsd$qldata[c(1,2)])
#' 
#' # augment design with two additional points
#' Xnew <- multiDimLHS(N=2,qsd$lower,qsd$upper,X=X,
#'            method="augmentLHS",type="matrix")
#' 
#' # new simulations
#' Xsim <- simQLdata(sim=qsd$simfn,nsim=10,X=Xnew)
#' 
#' # prepare data
#' Xdata <- setQLdata(Xsim,Xnew)
#' 
#' # do not re-estimate covariance parameters
#' qsd2 <- updateCovModels(qsd,Xdata,fit=FALSE) 
#'  
#' @seealso \code{\link{setQLdata}}, \code{\link{simQLdata}}, \code{\link{QLmodel}}
#' 
#' @author M. Baaske
#' @rdname updateCovModels
#' @export
updateCovModels <- function(qsd, nextData, fit = TRUE,
						cl = NULL, controls = list(), verbose=FALSE)
{		
	stopifnot(class(qsd) == "QLmodel")
	stopifnot(inherits(nextData,"QLdata"))	
	
	nnew <- NROW(nextData)	
	xdim <- attr(qsd$qldata,"xdim")
	nsim <- attr(qsd$qldata,"nsim")
	nsim.new <- attr(nextData,"nsim")
	
	nstat <- length(qsd$covT)	
	stid <- (xdim+1):(xdim+nstat)
	vid <- c((xdim+nstat+1):(xdim+2*nstat))
	
	vars <- qsd$qldata[vid]
	vars.new <- nextData[vid]
	
	# combine old and new data and
	# check columns also because L+ might not be given
	if(ncol(nextData) != ncol(qsd$qldata))
	  stop("The number of columns of argument `nextData` does not match with `qldata`.")
	
	# merge to new data, one (sample point) added
	qsd$qldata <- rbind(qsd$qldata,nextData)	
	Xs <- as.matrix(qsd$qldata[seq(xdim)])
	np <- nrow(Xs)
		
	if(length(controls)>0L) {		
		# set default optimization controls
		opts <- nloptr::nl.opts()
		opts[names(controls)] <- controls
	} else {
		# use stored optimization controls
		opts <- attr(qsd,"opts")	
	}		 
	# update function 
	fitit <- (fit && !(nrow(Xs) %% qsd$nfit))
	
	#c(xm$fix.nugget, xm$ptol*data[[i]][-(1:(np-nnew))] )
	.update <- function(covT, data, vars.new=NULL){
		mod <- lapply(1:length(covT),		
				function(i) {		   
					xm <- covT[[i]]
					# set starting point
					xm$start <- xm$param[xm$free]	
					if(!is.null(xm$fix.nugget)) {					  
					  xm$fix.nugget <-
						if(!is.null(vars.new)){							
						  c(xm$fix.nugget,vars.new[[i]])
						} else c(xm$fix.nugget,rep(xm$fix.nugget[1],nnew)) # re-use first 
		   			} # else (not using simulation variance for REML)
					if(fitit) {
						# store data for REML (remove afterwards in 'doREMLfit')
						xm$dataT <- data[[i]]										      			
					}
					xm
				})	
		if(fitit) {  				
			res <- doInParallel(mod, doREMLfit, Xs=Xs, opts=opts,
					 cl=cl, verbose=verbose)		
			if(!inherits(res,"error")) {
				.extractCovModels(res,verbose)
			} else {
				msg <- paste0("Error fitting covariance parameters.")
				message(msg)
				structure("message"=msg,"error"=res)				
			}	
		} else structure(mod, class = "krige")
	}
	
	tryCatch({
	  qsd$covT <- .update(qsd$covT,
			              qsd$qldata[stid],
						  nextData[vid]/nsim.new)
				  
	  # update kriging VARIANCE models
 	  # Cholesky terms are the data
	  if(qsd$var.type == "kriging"){
		if(is.null(qsd$covL))
		   stop("A covariance model for kriging the variance matrix must be defined but is `Null`.")
	    qsd$covL <-
	     if(attr(qsd$qldata,"Nb") > 0){ 
			 # if bootstrrapping nugget variances
			 # and only use simulations variances at new points
			 .update(qsd$covL,
					 qsd$qldata[grep("^L[^b]",names(qsd$qldata))],
					 nextData[grep("^Lb",names(qsd$qldata))])
		 } else {
			 .update(qsd$covL,qsd$qldata[grep("^L[^b]",names(qsd$qldata))],NULL) 
		 }	 
  	  }
 	}, error = function(e) {
	     msg <- .makeMessage("Failed to update covariance models: ",
				 conditionMessage(e))		
		return(.qleError(message=msg,call=match.call(),error=e))
	})

	return( qsd )
}
mbaaske/qle documentation built on Feb. 3, 2018, 11:02 a.m.