## ----Loading_library,include=FALSE--------------------------------------------
library(RGeostats)
constant.define("ASP",1)


## ----Create_the_Grid----------------------------------------------------------
ndim = 2
dx = c(1,1)
nx = c(300,250)
grid = db.create(nx=nx,dx=dx)


## ----Create_Model-------------------------------------------------------------
percentage = 20
extends = grid$extends
range   = extends[1] * percentage / 100
model   = model.create(vartype="K-Bessel",range=range, sill=1)


## ----Plot_Model---------------------------------------------------------------
plot(model,xlim=c(0,200),title="The Matern Model")


## ----Create_Neighborhood------------------------------------------------------
neigh = neigh.create(ndim=ndim,type=0)


## ----Non_conditional_simulations_using_1000_Turning_Bands---------------------
seed   = 432431
nbtuba = 1000
result.tb = simtub(,grid,model,seed=seed,nbtuba=nbtuba)
title = paste0("Non Conditional Simulation (",nbtuba," bands)")
plot(result.tb,pos.legend=1,title=title)


## ----Non_conditional_simulations_using_10000_Turning_Bands,eval=FALSE---------
## seed   = 432431
## nbtuba = 10000
## result.tb = simtub(,grid,model,seed=seed,nbtuba=nbtuba)
## title = paste0("Non Conditional Simulation (",nbtuba," bands)")
## plot(result.tb,pos.legend=1,title=title)


## ----Non_conditional_simulation_using_Regular_system--------------------------
seed   = 31415
nbsimu = 10
gext   = c(range,range)
result.nc = spde(,grid,model,seed=seed,gext=gext,nbsimu=nbsimu)
title = paste0("SPDE: one (out of ",nbsimu,") simulations")
plot(result.nc,pos.legend=1,title=title)


## ----Creating_irregular_data_set----------------------------------------------
ndata  = 200
data   = db.create(x1=extends[1]*runif(ndata),
                   x2=extends[2]*runif(ndata))
nbtuba = 1000
data   = simtub(,data,model,nbtuba=nbtuba)
data   = db.rename(data,"Simu.V1.S1","data")
title  = paste0("Data file (",ndata," samples)")
plot(data,name.post=1,pch=19,title=title)


## ----Standard_Kriging---------------------------------------------------------
result.krige = kriging(data,grid,model,neigh,uc="")


## ----Plotting_Standard_Kriging_Estimate---------------------------------------
title = paste0("Kriging Estimate") 
plot(result.krige,name="Kriging.data.estim",title=title)
plot(data,add=T,bg="white",pch=21)


## ----Plotting_Standard_Kriging_Stdev------------------------------------------
title = paste0("Kriging Standard Deviation")
plot(result.krige,name="Kriging.data.stdev",title=title)
plot(data,add=T,name.post=1,bg="white")


## ----Meshing_with_nqQ---------------------------------------------------------
title = paste0("Triangulation with switch 'nqQ'")
a = meshing(data,triswitch="nqQ",title=title,flag.point=FALSE,
            flag.plot=TRUE)
print(a)


## ----Meshing_with_nqQa50------------------------------------------------------
title = paste0("Triangulation with switch 'nqQa50'") 
a = meshing(data,triswitch="nqQa50",title=title,flag.point=FALSE,
            flag.plot=TRUE)
print(a)


## ----Meshing_with_extension---------------------------------------------------
title = paste0("Triangulation with extended area") 
a = meshing(data,triswitch="nqQ",dbaux=grid,gext=gext,
    title = title, flag.point=FALSE, flag.plot=TRUE, verbose=TRUE)
print(a)


## ----Performing_SPDE_Kriging_without_extension--------------------------------
result.kspde1 = spde(data,grid,model,flag.est=TRUE,flag.std=TRUE)


## ----Plotting_SPDE_Kriging_Estimate-------------------------------------------
title = paste0("Kriging Estimate (SPDE)")
plot(result.kspde1,name="SPDE.data.estim",title=title)
plot(data,add=TRUE,bg="white",pch=21)


## ----Plotting_SPDE_Kriging_Stdev----------------------------------------------
title = paste0("Kriging Standard Deviation (SPDE)")
plot(result.kspde1,name="SPDE.data.stdev",title=title)
plot(data,add=T,name.post=1,bg="white")


## ----Comparing_Kriging_Estimates_without_extension----------------------------
title = paste0("Comparing Estimations")
dummy = correlation(result.krige,"*estim","*estim",result.kspde1,title=title,
                    flag.iso=T,flag.same=T,flag.diag=T,
                    xlab="Traditional estimation",ylab="SPDE estimation")


## ----Comparing_Kriging_Stdev_without_extension--------------------------------
title = paste0("Comparing Standard Deviations")
dummy = correlation(result.krige,"*stdev","*stdev",result.kspde1,title=title,
                    flag.iso=T,flag.same=T,flag.diag=T,
                    xlab="Traditional estimation",ylab="SPDE estimation")


## ----Plotting_absolute_difference_of_Kriging_Estimates------------------------
result.kspde1 = db.add(result.kspde1,
diff.K = abs(result.kspde1[,"SPDE.data.estim"] -
             result.krige[,"Kriging.data.estim"]))
title = paste0("Difference of Estimations")
plot(result.kspde1,pos.legend=1,zlim=c(0,0.67),title=title)


## ----Plotting_absolute_difference_of_Kriging_Stdev----------------------------
result.kspde1 = db.add(result.kspde1,
diff.STD = abs(result.kspde1[,"SPDE.data.stdev"] -
               result.krige[,"Kriging.data.stdev"]))
title = paste0("Difference of Standard Deviations")
plot(result.kspde1,pos.legend=1,zlim=c(0,0.16),title=title)


## ----Performing_SPDE_Kriging_with_extension-----------------------------------
result.kspde2 = spde(data,grid,model,flag.est=TRUE,
                     flag.std=TRUE,gext=gext)


## ----Comparing_Kriging_Estimates_with_extension-------------------------------
title = paste0("Comparing Estimations")
dummy = correlation(result.krige,"*estim","*estim",result.kspde2,title=title,
                    flag.iso=T,flag.same=T,flag.diag=T,
                    xlab="Traditional estimation",ylab="SPDE estimation")


## ----Comparing_Kriging_Stdev_with_extension-----------------------------------
title = paste0("Comparing Standard Deviations")
dummy = correlation(result.krige,"*stdev","*stdev",result.kspde2,title=title,
                    flag.iso=T,flag.same=T,flag.diag=T,
                    xlab="Traditional estimation",ylab="SPDE estimation")


## ----Traditional_Conditional_simulation---------------------------------------
seed   = 543243
nbtuba = 1000
nbsimu = 1
result.tbc = simtub(data,grid,model,neigh,uc="",seed=seed,
                    nbsimu=nbsimu,nbtuba=nbtuba)


## ----Plotting_Conditional_simulation------------------------------------------
title = paste0("Conditional Simulation (",nbtuba," bands)",sep="")
plot(result.tbc,pos.legend=1,title=title)


## ----Performing_SPDE_conditional_simulations----------------------------------
nbsimu = 50
result.spc = spde(data,grid,model,seed=seed,gext=gext,
                  flag.modif=TRUE, nbsimu=nbsimu)


## ----Plotting_Mean_of_conditional_simulations---------------------------------
title = paste0("Average of ",nbsimu," simulations")
plot(result.spc,pos.legend=1,name="SPDE.mean",title=title)


## ----Plotting_Dispersion_of_conditional_simulations---------------------------
title = paste0("Standard Deviation of ",nbsimu," simulations")
plot(result.spc,pos.legend=1,name="SPDE.stdev",title=title)


## ----Comparing_Kriging_Estimate_and_Mean_of_conditional_simulations-----------
title = paste0("Comparing Kriging and Mean of Simulations")
dummy = correlation(result.krige,"*estim","*mean",result.spc,title=title,
                    flag.iso=T,flag.same=T,flag.diag=T,
                    xlab="Direct Estimation",ylab="From Simulations")


## ----Comparing_Kriging_Stdev_and_Dispersion_of_conditional_simulations--------
title = paste0("Comparing Standard Deviation and Dispersion of Simulations")
dummy = correlation(result.krige,"*stdev","*stdev",result.spc,title=title,
                    flag.iso=T,flag.same=T,flag.diag=T,
                    xlab="Direct Estimation",ylab="From Simulations")


## ----Creating_Large_data_set--------------------------------------------------
ndata.big = 10000
data.big  = db.create(x1=extends[1]*runif(ndata.big),
                      x2=extends[2]*runif(ndata.big))
nbtuba    = 1000
data.big  = simtub(,data.big,model,nbtuba=nbtuba)
data.big  = db.rename(data.big,"Simu.V1.S1","data")


## ----Plotting_Large_data_set--------------------------------------------------
title = paste0("Large data set (",ndata.big," samples)")
plot(data.big,name.post=1,cex=0.15,title=title)


## ----Kriging_Large_data_set---------------------------------------------------
result.kspde.big = spde(data.big,grid,model,flag.est=TRUE,flag.std=TRUE)


## ----Plotting_Kriging_Estimate_with_Large_data_set----------------------------
title = paste0("Estimation (SPDE) with large data set") 
plot(result.kspde.big,name="SPDE.data.estim",title=title)
plot(data.big,add=T,name.post=1,cex=0.1)


## ----Plotting_Kriging_Stdev_with_Large_data_set-------------------------------
title = paste0("Standard Deviation (SPDE) with large data set")
plot(result.kspde.big,name="SPDE.data.stdev",title=title)
plot(data.big,add=T,name.post=1,cex=0.1)


## ----Performing_Conditional_Simulations_with_Large_data_set-------------------
nbsimu = 10
result.spc.big = spde(data.big,grid,model,seed=seed,gext=gext,
                      nbsimu=nbsimu)
title = paste0("Conditional simulation (one out of ",nbsimu,") with large data set")
plot(result.spc.big,pos.legend=1,title=title)