Nothing
# This function estimates any of the the NA parameters, or returns c and d if all parameters are given
CKrigcoefs <- function(xy,
nuniquerows, ## as given by selectFn
covfnparamA=NA, # excluding smoohtness, but output covfnparam includes it.
lambdaA=NA,
minSmoothness=blackbox.getOption("minSmoothness"),
maxSmoothness=blackbox.getOption("maxSmoothness"),
miscOptions=blackbox.getOption("miscOptions"),
initCovFnParam=NULL,
verbosity=blackbox.getOption("verbosity"),
optimizers=blackbox.getOption("optimizers")
){
nrowxy <- nrow(xy)
ncolxy <- ncol(xy)
# Test on length(covfnparamA) assumes that when scale params are known, smoothness is too
# (=> testing smoothness range is irrelevant).
# In standard usage either covfnparamA is zero-length (and lambda is unknown)
# => the estimation of cov params is performed
# or all parameters (including lambda) are known and no optimization is performed.
# It is also possible that the user sets a scale parameter, and and smoothness through min/maxSmoothness,
# but lambda still needs to be estimated: this is a subcase of the optimize case.
# (we might only test lambda to distinguish the cases)
# The user is currently not allowed to set all correlation params through the CovFnParam vector:
# input CovFnParam should contain scale parameters only.
# The user is also not allowed to provide lambda (through blackbox.getOption("lambdaEst")).
varCorrpars <- (length(covfnparamA)< ncolxy ## shorter than scale params+smoothness
|| any(is.na(covfnparamA)))
optimise <- (varCorrpars || is.na(lambdaA))
method <- intersect(optimizers,c("bobyqa","L-BFGS-B","lbfgsb3c")) ## non-default methods
if(length(method)==0L) method <- "NLOPT_LN_BOBYQA" ## default for this part of code.
if(length(method)!=1L) stop("length(method)!=1L in CKrigcoefs") ## incompatible with switch below
verbosityFromObjective <- as.integer((optimise && "L-BFGS-B" %in% method))*verbosity ## FR->FR should control in the C++ code using the eval counter
verbosityFromoptimizer <- switch(method,
"bobyqa"=round(2*as.integer(optimise)*verbosity), ## default=2
"L-BFGS-B"=0, ## verbosity from optim() is awkward
as.integer(optimise)*verbosity
)
if ("Migraine" %in% blackbox.getOption("usedBy")) {
init_factor <- 2/3
} else init_factor <- 2/30 ## better when rugged landscapes are possible
success <- newCSmooth(xy= t(xy), ## newCSmooth is an Rcpp export
nrowxy=nrowxy,
ncolxy=ncolxy,
nuniquerows=nuniquerows,
GCV=0, # 0 for GCV, 1 for match_fs2hat_pure_error
optimiseBool=optimise,
verbosity=verbosityFromObjective)
if ( ! success ) {
resu <- list()
## go directly to deleteCSmooth
} else if ( ! optimise ) {
resu <- Krig_coef_Wrapper(covfnparamA,lambdaA)[c("c","d","CKrigidx")] # throwing away u and D
# do not deleteCSmooth(), CSmooth pointer stored in CKrigptrTable
} else { # optimize
if (varCorrpars) { ## need to estim scale params and/or smoohtness (not sure that it can estimate only some of them)
if ( is.null(initCovFnParam) ) initCovFnParam <- rep(0, ncolxy)
xx <- xy[,-ncolxy,drop=FALSE]
KgLow <- sapply(xx, min)
KgUp <- sapply(xx, max)
maxrange <- KgUp-KgLow
GCVlowerFactor <- 20 # or 20 ? cf comments on Migraine default in Migraine code
GCVupperFactor <- 5
lower <- maxrange/(GCVlowerFactor*nuniquerows)
upper <- maxrange*GCVupperFactor
init_factor <- rep(init_factor,length(lower))
minSmoothness <- min(maxSmoothness,max(minSmoothness,1.001)); # >1, cf def Matern
varSmoothness <- (maxSmoothness-minSmoothness>1e-6)
if (varSmoothness) { # Variable Smoothness case
maxrange <- c(maxrange,maxSmoothness-minSmoothness)
lower <- c(lower,minSmoothness)
upper <- c(upper,maxSmoothness)
init_factor <- c(init_factor,2/3)
fixedSmoothness <- numeric(0) ## not c() which is NULL...
} else {
initCovFnParam <- initCovFnParam[seq_len(length(lower))]
fixedSmoothness <- maxSmoothness
}
for (it in seq_len(length(lower))) {
if (initCovFnParam[it]<lower[it] ## includes default case
|| initCovFnParam[it]>upper[it]) {
initCovFnParam[it]=lower[it]+maxrange[it]*init_factor[it]
}
}
if ("L-BFGS-B" %in% method) {
control <- list(parscale=(upper-lower), trace=verbosityFromoptimizer)
optr <- optim(par=initCovFnParam,fn=GCV_lamVar_covFix_Wrapper,method="L-BFGS-B",
lower=lower,upper=upper,control=control,
fixedSmoothness=fixedSmoothness,returnFnvalue=TRUE)
solution <- optr$par
} else if ("lbfgsb3c" %in% method){ ## very slow
control <- list(trace=verbosityFromoptimizer)
if ( ! requireNamespace("lbfgsb3c",quietly=TRUE) ) {
stop("Package lbfgsb3c not installed.")
}
optr <- lbfgsb3c::lbfgsb3c(par=initCovFnParam,fn=GCV_lamVar_covFix_Wrapper,lower=lower,upper=upper,
control=control,
fixedSmoothness=fixedSmoothness,returnFnvalue=TRUE)
solution <- optr$par # change for lbfgsb3c
} else if ("bobyqa" %in% method){ ## marginally better than optim ?
if ( ! requireNamespace("minqa",quietly=TRUE) ) {
stop("Package minqa not installed.")
}
control <- list(rhobeg=min(abs(upper-lower))/20,iprint=verbosityFromoptimizer)
control$rhoend <- max(1,control$rhobeg)/1e6
optr <- minqa::bobyqa(par=initCovFnParam,fn=GCV_lamVar_covFix_Wrapper,lower=lower,upper=upper,
control=control,
fixedSmoothness=fixedSmoothness,returnFnvalue=TRUE)
solution <- optr$par
} else if ("NLOPT_LN_BOBYQA" %in% method){ ## may be the fastest
optr <- nloptr(x0=unlist(initCovFnParam),eval_f=GCV_lamVar_covFix_Wrapper,lb=lower,ub=upper,
opts=list(algorithm="NLOPT_LN_BOBYQA",xtol_rel=1.0e-4,maxeval=-1,print_level=verbosityFromoptimizer),
fixedSmoothness=fixedSmoothness,returnFnvalue=TRUE)
solution <- optr$solution
} else stop("Unknown 'method' in CKrigcoefs()")
solution <- c(solution,fixedSmoothness) # required when minSmoothness=maxSmoothness
resu <- list(covfnparam=solution, # full including smoothness whether fixed or estimated
fnEvalCount=getFnEvalCount(),method=method)
} else { ## no estim of scale nor smoohtness but we will estimate lambda
resu <- list(covfnparam=c(covfnparamA,maxSmoothness)) # also including smoothness
}
if (is.na(lambdaA)) resu$lambda <- GCV_lamVar_covFix_Wrapper(resu$covfnparam,fixedSmoothness=numeric(0),returnFnvalue=FALSE)
deleteCSmooth() ## this code conditional on optimisation
}
return(resu) # is a list with elements depending on call:
# if optimise covfnparam -> returns covfnparam, lambda, fnEvalCount, method
# else if optimise lambda -> returns (input covfnparamA), lambda
# else returns c("c","d","CKrigidx")
} ## end def CKrigcoefs
# R-style implementation of default values
R_GCV_lamVar_covFix <- function(a,fixedSmoothness=numeric(0L),returnFnvalue=TRUE) {
## call to R-to-C wrapper with explicit values
GCV_lamVar_covFix_Wrapper(a=a,fixedSmoothness=fixedSmoothness,returnFnvalue=returnFnvalue)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.