R/scalingFun.R

In IRSN/DiceKriging: Kriging Methods for Computer Experiments

##===============================================================
## 'scalingFun1d' computes the 'scaling' transformation at
## the one-dimensional design points x.
##
## The knots and (positive) derivatives definitng the transfo-
## -mation are given in numeric vectors 'knots' and 'eta'.
##
##===============================================================

#' @example plot(function(x)scalingFun1d(x,c(0,1),c(.5,.5)))
#' @example plot(function(x)scalingFun1d(x,c(0,1),c(.5,.75)))
#' @example plot(function(x)scalingFun1d(x,0,.5))
#' @example plot(function(x)scalingFun1d(x,c(0,1),c(.5,.75)),xlim=c(-1,2))
#' @example plot(function(x)scalingFun1d(x,c(0,1.5),c(.5,1.75)),xlim=c(-1,2))
#' @example plot(function(x)scalingFun1d(x,c(0,.5,1),c(10.5,.5,1.75)),xlim=c(-1,2))
scalingFun1d <- function(x, knots, eta){

    n <- length(x)
    nKnots <- length(knots)

    #if ( any(x < knots[1]-1E-6) | any(x > knots[nKnots]+1E-6) )
    #    warning("'x' values should be inside the knots (otherwise using closest knot)\nknots =",paste(knots,collapse=","),"\nx = ",paste(x,collapse=","))

    ix_upper= which(x > knots[nKnots])
    ix_lower = which(x < knots[1])
    if (length(ix_lower)>0 & length(ix_upper)>0)
        ix_inside = (1:n)[-c(ix_upper,ix_lower)]
    else if (length(ix_lower)>0 & length(ix_upper)==0)
        ix_inside = (1:n)[-ix_lower]
    else if (length(ix_upper)>0 & length(ix_lower)==0)
        ix_inside = (1:n)[-ix_upper]
    else ix_inside=1:n

    newscale_inside = scalingFun1d.inside(x[ix_inside],knots, eta)

    newscale <- rep(NA, n)
    newscale[ix_inside] = newscale_inside

    ## Support for outside knots
    if(length(ix_lower)>0) newscale[ix_lower]=eta[1]*(x[ix_lower]-knots[1])

    if(length(ix_upper)>0) {
        y0=scalingFun1d.inside(knots[nKnots],knots, eta)
        newscale[ix_upper]=eta[nKnots]*(x[ix_upper]-knots[nKnots])+y0
    }

    return(newscale)
}


#' Do NOT support when x is outside bounds of knots...
#' @perf x=runif(100); k=seq(0,1,,10);e=runif(10); scalingFun1d.inside.faster(x,k,e)
#' @perf Rprof("scaling.prof", line.profiling=TRUE); x=runif(100); k=seq(0,1,,10);e=runif(10); for(i in 1:10000) scalingFun1d.inside(x,k,e); Rprof(NULL); summaryRprof("scaling.prof", lines = "show")
scalingFun1d.inside <- function(x, knots, eta){

    n <- length(x)
    nKnots <- length(knots)

    if (nKnots == 1) return(eta*(x-knots))

    if ( any(x < knots[1]) | any(x > knots[nKnots]) )
        stop("'x' values should be inside the knots (otherwise using closest knot)\nknots =",paste(knots,collapse=","),"\nx = ",paste(x,collapse=","))

    if (length(x)>1){
        xs <- sort(x, index.return = TRUE)
        xsorted <- xs$x
        ind <- xs$ix
    } else{
        xsorted=x
        ind=1
    }

    scale <- rep(0, n)

    ## cat("Calling C\n")
    res <- .C("Scale",
              n = as.integer(n),
              nKnots = as.integer(nKnots),
              x = as.double(xsorted),
              knots = as.double(knots),
              eta = as.double(eta),
              scale = as.double(scale))

    ## CAUTION here: the values are for SORTED x values
    ## they must be put back in the original order.
    newscale <- rep(NA, n)
    newscale[ind] <- res$scale

    return(newscale)
}

##==================================================================
## 'scalingFun' applies d one-dimensional 'scaling' transformations
## to n d-dimensional design points
##
## o 'X' must be a n*d matrix,
##
## o 'knots' and 'eta' are list of lentgh d. In both cases,
##    the element i contains the knots/derivatives for the
##    dimension i
##
## At this time, no control on the list/matrix dimension or length
##
##==================================================================

scalingFun <- function(X, knots, eta, plot = FALSE) {

    d <- NCOL(X)

    transX <- matrix(NA, nrow = NROW(X), ncol = NCOL(X))
    dimnames(transX) <- dimnames(X)
    X <- as.matrix(X)

    for (i in 1:d){
        name = i; if(!is.null(colnames(X))) name = colnames(X)[i]
        transX[ , i] <- scalingFun1d(x = X[ , i], knots = knots[[name]], eta = eta[[name]])
    }

    if (plot) {
        for(i in 1:d) {
            plot(X[ ,i], transX[ ,i], type = "o",
                 pch = 21, col = "orangered", cex = 0.8,
                 xlab = "x", ylab = "x 'scaled'")
            abline(v = knots[[i]], col = "SpringGreen3")
        }

    }

    return(transX)

}


affineScalingFun <- function(X, knots, eta) {
    # Here X is meant to be a n*d matrix,
    # knots is a (K+1) vector (in this special case),
    # and eta is a d*(K+1) matrix of coefficients.

    d <- dim(X)[2]

    transformed_X <- matrix(NA, nrow=nrow(X), ncol=ncol(X))

    for(i in seq(1,d)){
        transformed_X[,i] <- scalingFun1d(x=X[,i], knots=knots, eta = eta[i,])
    }

    return(transformed_X)

}