# jn_optim_parallel2: Parallel jn criterion In KrigInv: Kriging-Based Inversion for Deterministic and Noisy Computer Experiments

## Description

Evaluation of the parallel jn criterion for some candidate points, assuming that some other points are also going to be evaluated. To be used in optimization routines, like in `max_sur_parallel`. To avoid numerical instabilities, the new points are evaluated only if they are not too close to an existing observation, or if there is some observation noise. The criterion is the integral of the posterior expected `"jn"` uncertainty, which is the posterior expected variance of the excursion set's volume.

## Usage

 ```1 2 3 4 5``` ```jn_optim_parallel2(x, other.points, integration.points, integration.weights = NULL, intpoints.oldmean, intpoints.oldsd, precalc.data, model, T, new.noise.var = NULL, batchsize, current.sur,ai_precalc) ```

## Arguments

 `x` Input vector of size d at which one wants to evaluate the criterion. This argument corresponds to only ONE point. `other.points` Vector giving the other `batchsize-1` points at which one wants to evaluate the criterion `integration.points` Matrix of points for numerical integration. Two cases are handled. If the `"jn"` importance sampling distribution has been used, this is a (2p)*d matrix containing p couples of integration points in D x D, where D is the integration domain. If instead a p*d matrix is given, then the function will perform the integration by using all `p^2` couples. `integration.weights` Vector of size p corresponding to the weights of these integration points. `intpoints.oldmean` Vector of size p, or 2p, corresponding to the kriging mean at the integration points before adding `x` to the design of experiments. `intpoints.oldsd` Vector of size p, or 2p, corresponding to the kriging standard deviation at the integration points before adding `x` to the design of experiments. `precalc.data` List containing useful data to compute quickly the updated kriging variance. This list can be generated using the `precomputeUpdateData` function. `model` Object of class `km` (Kriging model). `T` Array containing one or several thresholds. `new.noise.var` Optional scalar value of the noise variance of the new observations. `batchsize` Number of points to sample simultaneously. The sampling criterion will return batchsize points at a time for sampling. `current.sur` Current value of the sur criterion (before adding new observations) `ai_precalc` When multiple thresholds are used (i.e. when T is a vector), this is an nT*p matrix with ith row equal to `intpoints.oldmean-T[i]`. The argument DOES ALWAYS need to be filled, even if there is only one threshold T.

## Details

The first argument `x` has been chosen to be a vector of size d so that an optimizer like genoud can optimize it easily. The second argument `other.points` is a vector of size (batchsize-1)*d corresponding to the batchsize-1 other points.

## Value

Parallel jn value

## Author(s)

Clement Chevalier (University of Neuchatel, Switzerland)

## References

Chevalier C., Bect J., Ginsbourger D., Vazquez E., Picheny V., Richet Y. (2014), Fast parallel kriging-based stepwise uncertainty reduction with application to the identification of an excursion set, Technometrics, vol. 56(4), pp 455-465

Chevalier C., Ginsbourger D. (2014), Corrected Kriging update formulae for batch-sequential data assimilation, in Pardo-Iguzquiza, E., et al. (Eds.) Mathematics of Planet Earth, pp 119-122

`EGIparallel`, `max_sur_parallel`
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73``` ```#jn_optim_parallel2 set.seed(9) N <- 20 #number of observations T <- c(80,100) #threshold or thresholds testfun <- branin #a 20 points initial design design <- data.frame( matrix(runif(2*N),ncol=2) ) response <- testfun(design) #km object with matern3_2 covariance #params estimated by ML from the observations model <- km(formula=~., design = design, response = response,covtype="matern3_2") ###we need to compute some additional arguments: #integration points, and current kriging means and variances at these points integcontrol <- list(n.points=200,distrib="jn",init.distrib="MC") obj <- integration_design(integcontrol=integcontrol,lower=c(0,0),upper=c(1,1), model=model,T=T) integration.points <- obj\$integration.points integration.weights <- obj\$integration.weights pred <- predict_nobias_km(object=model,newdata=integration.points, type="UK",se.compute=TRUE) intpoints.oldmean <- pred\$mean ; intpoints.oldsd<-pred\$sd #another precomputation precalc.data <- precomputeUpdateData(model,integration.points) nT <- 2 # number of thresholds ai_precalc <- matrix(rep(intpoints.oldmean,times=nT), nrow=nT,ncol=length(intpoints.oldmean),byrow=TRUE) ai_precalc <- ai_precalc - T # substracts Ti to the ith row of ai_precalc batchsize <- 4 other.points <- c(0.7,0.5,0.5,0.9,0.9,0.8) x <- c(0.1,0.2) #one evaluation of the jn_optim_parallel criterion2 #we calculate the expectation of the future "sur" uncertainty when #1+3 points are added to the doe #the 1+3 points are (0.1,0.2) and (0.7,0.5), (0.5,0.9), (0.9,0.8) jn_optim_parallel2(x=x,other.points,integration.points=integration.points, integration.weights=integration.weights, intpoints.oldmean=intpoints.oldmean,intpoints.oldsd=intpoints.oldsd, precalc.data=precalc.data,T=T,model=model, batchsize=batchsize,current.sur=Inf,ai_precalc=ai_precalc) n.grid <- 20 #you can run it with 100 x.grid <- y.grid <- seq(0,1,length=n.grid) x <- expand.grid(x.grid, y.grid) jn_parallel.grid <- apply(X=x,FUN=jn_optim_parallel2,MARGIN=1,other.points, integration.points=integration.points, integration.weights=integration.weights, intpoints.oldmean=intpoints.oldmean,intpoints.oldsd=intpoints.oldsd, precalc.data=precalc.data,T=T,model=model, batchsize=batchsize,current.sur=Inf,ai_precalc=ai_precalc) z.grid <- matrix(jn_parallel.grid, n.grid, n.grid) #plots: contour of the criterion, doe points and new point image(x=x.grid,y=y.grid,z=z.grid,col=grey.colors(10)) contour(x=x.grid,y=y.grid,z=z.grid,15,add=TRUE) points(design, col="black", pch=17, lwd=4,cex=2) points(matrix(other.points,ncol=2,byrow=TRUE), col="red", pch=17, lwd=4,cex=2) i.best <- which.min(jn_parallel.grid) points(x[i.best,], col="blue", pch=17, lwd=4,cex=3) #plots the real (unknown in practice) curve f(x)=T testfun.grid <- apply(x,1,testfun) z.grid.2 <- matrix(testfun.grid, n.grid, n.grid) contour(x.grid,y.grid,z.grid.2,levels=T,col="blue",add=TRUE,lwd=5) title("Contour lines of jn_parallel criterion (black) and of f(x)=T (blue)") ```