## ----Loading_library,include=FALSE--------------------------------------------
rm(list = ls())
library(RGeostats)
library(resample)
knitr::opts_chunk$set(echo = TRUE, fig.align="center")
constant.define("asp",1)
set.seed(32432)


## ----Creating the Model and the Variables on Grid-----------------------------
A       <- matrix(c(1,0.8,0.2,0.8,1,0.7,0.2,0.7,1),ncol=3,nrow=3,byrow=T)
a.sill  <- t(A) %*% A
a.model <- model.create(vartype=3,nvar=3,sill=a.sill,range=10)
a.model <- model.create(vartype=1,nvar=3,model=a.model,sill=diag(1:3)/2)

# Simulation du modele
nbtuba = 500
a.grid <- db.create(flag.grid=T,x0=c(0,0),nx=c(100,100),dx=c(1,1))
a.grid <- simtub(dbout=a.grid,model=a.model,nbtuba = nbtuba,seed=4143)


## ----Comparing empirical variograms and the Model-----------------------------
plot(vario.grid(a.grid),title="Model of Raw Variables")
plot(a.model, add=T)
dev.off()


## ----Model reproduction-------------------------------------------------------
round(db.stat(a.grid, fun = "var", names = 4:6, flag.mono = F),2)
a.sill


## ----Computing PCA with RGeostats---------------------------------------------
a.pca <- pca.calc(a.grid)


## ----Extracting Data----------------------------------------------------------
Z <- as.matrix(a.grid[,"Simu*"])


## ----Mean comparison----------------------------------------------------------
# Explicit calculation
means = colMeans(Z)
means
# Contents of the PCA object
a.pca$mean


## ----Standard deviation comparison--------------------------------------------
# Explicit calculation
sigmas = colStdevs(Z)
sigmas
# Contents of the PCA object
a.pca$sigma
N = dim(Z)[1]
sqrt(a.pca$sigma * a.pca$sigma * N / (N-1)) 


## ----Eigen values comparison--------------------------------------------------
# Explicit calculation
Z_s = scale(Z)
Eig = eigen(var(Z_s))
Eig
# Contents of the PCA object
a.pca$eigen
a.pca$pcaz2f


## ----Transforming Variables into Factors using RGeostats----------------------
a.grid_1 <- pca.z2f(a.grid,a.pca)
Y_1 = as.matrix(a.grid_1[,"Factor*"])


## ----Transforming Variables into Factors using Explicit calculations----------
P0      <- Eig$vectors
lambda0 <- Eig$values
# Building factors
Y_2 <- Z_s %*% P0


## ----Comparing Factors--------------------------------------------------------
a.sign = P0 / a.pca$pcaz2f
a.sign = apply(a.sign,2,mean)
a.sign
Y_1 = sweep(Y_1,2,a.sign,FUN="*")
apply(Y_1 - Y_2,2,range)


## ----Graphic comparison-------------------------------------------------------
a.grid_2 <- pca.f2z(a.grid_1,a.pca)
res.0 <- a.grid_2[,"Simu*"]
res.1 <- a.grid_2[,"Variable*"]
res.2 <- Y_2 %*% t(P0)
res.2 = sweep(res.2,2,sigmas,FUN="*")
res.2 = sweep(res.2,2,means,FUN="+")
plot(res.2[,1],res.1[,1],col=1,pch=20,cex=0.1,
     xlab ="input data (Z)", 
     ylab="Back-transformation of the factors",xlim=c(-6,6),asp=1,
     main="Comparison between Initial Variables and back-transformed Factors")
points(res.1[,1],res.0[,1],col=2,pch=20,cex=0.1)
abline(0,1)


## ----Using db.add.CL----------------------------------------------------------
#standardisation of the input variables
a.grid.std = db.normalize(a.grid,names="Simu*",oper="scal")
# Computing the factors with db.add.CL
a.grid_3 <- db.add.CL(a.grid.std,
                      Y.names = "Norm*", Z.names = paste("Y",1:3,sep="."),
                      mat.CL=t(a.pca$pcaz2f))
summary(a.grid_3$Y.1 - a.grid_1$Factor.1)
summary(a.grid_3$Y.2 - a.grid_1$Factor.2)
summary(a.grid_3$Y.3 - a.grid_1$Factor.3)


## ----Computing MAF using RGeostats--------------------------------------------
a.maf <- maf.calc(a.grid,h0=1,dh=0.01,dirvect = 0,tolang=0)


## ----Computing MAF using Explicit calculations--------------------------------
# Standardizing factors
D_inv_sqrt_lambda0 <- diag(1/sqrt(lambda0))
Y_2_s <- Y_2 %*% D_inv_sqrt_lambda0
# Spatial increments 
tmp  <- db.add(a.grid,as.data.frame(Y_2_s))
tmp0 <- db.reduce(db.sel(tmp,x1<99),names="V*",flag.keep.grid = T)
tmp1 <- db.reduce(db.sel(tmp,x1>0),names="V*",flag.keep.grid = T)
tmp0 <- as.matrix(db.extract(tmp0,"V*"))
tmp1 <- as.matrix(db.extract(tmp1,"V*"))
dY_sh <- tmp1-tmp0
# ACP of the Variance-Covariance of the increments
Eig <- eigen(var(dY_sh))
Ph  <- Eig$vectors
lambdah <- Eig$values
# Ordering MAFs by increasing eigen values
tmp <- order(lambdah)
Ph  <- Ph[,tmp]
lambdah <- lambdah[tmp]


## ----Comparing PCA Eigen values-----------------------------------------------
# Using maf calculations
a.maf$eigen

# Using explicit calculations
lambdah


## ----Comparing MAF Eigen values-----------------------------------------------
# Using maf calculations
a.maf$pcaz2f

# Using explicit calculations
P <- P0 %*% D_inv_sqrt_lambda0 %*% Ph
P


## ----Performing transformation using RGeostats--------------------------------
a.grid.xsi_1 <- pca.z2f(a.grid,a.maf,radix="MAF")
xsi_1 = db.extract(a.grid.xsi_1,"MAF*")


## ----Performing transformation using Explicit calculations--------------------
xsi_2 <- Z_s %*% P
a.sign = P / a.maf$pcaz2f
a.sign = round(apply(a.sign,2,mean),0)
a.sign
xsi_2 = sweep(xsi_2,2,a.sign,FUN="*")
a.grid.xsi_2 = db.add(a.grid, as.data.frame(xsi_2))


## ----Comparing factors--------------------------------------------------------
apply(xsi_1 - xsi_2,2,range)


## ----Maf models---------------------------------------------------------------
plot(vario.grid(a.grid.xsi_1),title="Model for MAF factors")
plot(vario.grid(a.grid.xsi_2),add=T)


## ----MAF representation-------------------------------------------------------
plot(a.grid.xsi_1,name="MAF.1")
plot(a.grid.xsi_1,name="MAF.2")
plot(a.grid.xsi_1,name="MAF.3")