## ----Loading_library,include=FALSE--------------------------------------------
library("RGeostats")
library(knitr)
rm(list=ls())
set.seed(13241)


## ----creating-data-nocorr, echo=FALSE, results='hide'-------------------------
# les points de données
data.cl  <- t(c(2, 4))
colnames(data.cl) <- c("V1", "V2")
rownames(data.cl) <- c("Z")
data.nbr <- 500
data.dx  <- c(10, 10)
data <- db.create(flag.grid = FALSE,
                  x1 = runif(data.nbr, min = 0, max = data.dx[1]),
                  x2 = runif(data.nbr, min = 0, max = data.dx[2]))
sel_data   <- sample(x = c(TRUE, FALSE), size = data.nbr,  
                     replace = TRUE, prob = c(2/3, 1/3))
sel_target <- (!sel_data)
data <- db.add(data, sel_data, sel_target)
data <- db.locerase(data, loctype = "z")
# données homotopiques sans corrélation
model.V <- model.join(
  model.create(vartype = 2, range = data.dx[1]/4), 
  model.create(vartype = 5, range = data.dx[1]/2)
)
data <- simtub(dbout = data, 
               model = model.CL(model.V, 
                                mat.CL = matrix(c(1,0),nrow=1,ncol=2)), 
               seed = 1244,
               nbsimu = 1, nbtuba = 1000)
data <- db.rename(data, name = "Simu.V1.S1", newname = "V1")
data <- simtub(dbout = data, 
               model = model.CL(model.V, 
                                mat.CL = matrix(c(0,1),nrow=1,ncol=2)), 
               seed = 1244/2,
               nbsimu = 1, nbtuba = 1000)
data <- db.rename(data, name = "Simu.V1.S1", newname = "V2")
data <- db.add.CL(data, Y.names = colnames(data.cl), 
                  Z.names = rownames(data.cl), mat.CL = data.cl)
data <- db.locerase(data, loctype = "z")
# données hétérotopiques
for (i in 1:2) {
  W_obs <- db.extract(data, names = paste0("V", i), flag.compress = FALSE)
  W_obs[sample(x = c(TRUE, FALSE), size = data$nech,  replace = TRUE,
               prob = c(1/2, 1/2))] <- NA
  data <- db.add(data, W_obs)
  data <- db.rename(data, name = "W_obs", 
                    newname = paste0("V", i, "_obs"))
}
data <- db.add.CL(data, Y.names = paste0("V", 1:2, "_obs"), 
                  Z.names = "Z_obs", mat.CL = data.cl)
data <- db.locerase(data, loctype = "z")
rm(W_obs)

# graphique des données utilisées dans le test
plot(db.locate(data, names = "sel_data", loctype = "sel", 
               flag.locnew = TRUE), col = "blue", pch = 1,
     title= "Données")
plot(db.locate(data, names = "sel_target", loctype = "sel", 
               flag.locnew = TRUE), col = "red", pch = 19, add = TRUE)
legend("bottomright", 
       c("data", "target"),
       col = c("blue", "red"),
       pch = c(1, 19))


## ----voisinage----------------------------------------------------------------
neigh <- neigh.create(ndim = 2, type = 0)


## ----def-function, echo=FALSE-------------------------------------------------
# Définition d'une fonction pour les différents cas considérés
QC_LC_V12 <- function(cas) {
  Z.nm <- paste0("Z", cas)
  rownames(data.cl) <- Z.nm
  #cokrigeage des valeurs V1, V2 à partir de V1 ou V2
  data <- db.locate(data, names = "sel_data", loctype = "sel", flag.locnew = TRUE)
  data <- kriging(
    dbin  = db.locate(data, names = paste0("V", 1:2, cas), loctype = "z", flag.locnew = TRUE),
    dbout = db.locate(data, names = "sel_target", loctype = "sel", flag.locnew = TRUE),
    model = model.V,
    neigh = neigh,
    calcul = "point", uc = "1", flag.est = TRUE, flag.std = TRUE, 
    radix = "ck"
  )
  
  data <- db.add.CL(data, 
                    Y.names = paste0("ck.V", 1:2, cas, ".estim"), 
                    Z.name  = paste0("ck2.", Z.nm, ".estim"), 
                    mat.CL = data.cl)
  data <- db.delete(data, names = paste0("ck.V*", cas, ".*"))
  
  #cokrigeage de la CL Z
  data <- db.locate(data, names = "sel_data", loctype = "sel", 
                    flag.locnew = TRUE)
  
  data <- kriging(
    dbin  = db.locate(data, names = paste0("V", 1:2, cas), 
                      loctype = "z", flag.locnew = TRUE),
    dbout = db.locate(data, names = "sel_target", loctype = "sel", 
                      flag.locnew = TRUE),
    model = model.V, mat.CL = data.cl,
    neigh = neigh,
    calcul = "point", uc = "1", flag.est = TRUE, flag.std = TRUE, 
    radix = "ck"
  )
  
  # krigeage de Z il faut que X1 et X2 soient définies
  data <- db.locate(data, names = "sel_data", loctype = "sel", 
                    flag.locnew = TRUE)
  data <- kriging(
    dbin  = db.locate(data, names = Z.nm, loctype = "z", 
                      flag.locnew = TRUE),
    dbout = db.locate(data, names = "sel_target", loctype = "sel",
                      flag.locnew = TRUE),
    model = model.CL(model.V, mat.CL = data.cl),
    neigh = neigh,
    calcul = "point", uc = "1", flag.est = TRUE, flag.std = TRUE, 
    radix = "ok"
  )
  
  # sélection des données de contrôles
  data <- db.locate(data, names = "sel_target", loctype = "sel", flag.locnew = TRUE)
  # comparaison CL(est) vs. est(CL)
  correlation(db=data, size0=0.05,
              paste("ck", Z.nm, "estim", sep = "."),
              paste("ck2", Z.nm, "estim", sep = "."),
              flag.aspoint=TRUE, flag.diag=TRUE, 
              flag.ce=TRUE, flag.iso=TRUE,
              diag.col="red", diag.lwd=2, 
              ce.type="b", ce.col="blue", ce.lty=2, ce.pch=19,
              xlab="Cokrigeage de la C.L.",
              ylab="C.L. des Krigeages",
              title="Estimations")
  
  # Comparaison des écarts-types de (co)-krigeage
  kable(
    db.stat.multi(data, 
                  names = paste(c("ck", "ok"), Z.nm, "stdev", sep=".")),
    caption = "Statistiques des écarts-types d'estimation")
  correlation(db=data, size0=0.05,
              paste("ck", Z.nm, "stdev", sep = "."), 
              paste("ok", Z.nm, "stdev", sep = "."), 
              flag.aspoint=TRUE, flag.diag=TRUE, 
              flag.ce=TRUE, flag.iso=TRUE,
              diag.col="red", diag.lwd=2, 
              ce.type="b", ce.col="blue", ce.lty=2, ce.pch=19,
              xlab="Cokrigeage de la C.L.",
              ylab="C.L. des Krigeages",
              title="Ecart-types d'estimation")
  
  # comparaison Z vs. Z*
  est.list <- c("ck", "ok")
  var.nm   <- "Z"
  for (est.item in est.list) {
    if (est.item == "ck") {
      title = paste("Nuage de corrélation Z|Z* - Cokrigeage")
    } else {
      title = paste("Nuage de corrélation Z|Z* - Krigeage")
    }
    est.nm = paste(est.item, Z.nm, "estim", sep = ".")
    correlation(db=data, size0=0.05,
                est.nm, var.nm,
                flag.aspoint=TRUE, flag.diag=TRUE, 
                flag.ce=TRUE, flag.iso=TRUE,
                diag.col="red", diag.lwd=2, 
                ce.type="b", ce.col="blue", ce.lty=2, ce.pch=19,
                xlab="Estimation", ylab="Réalité", title=title)
  }
  data
}

# Calcul des scores de validation Z vs. Z* (cf. Thèse de F. FOUEDJIO p. 177)
db.score <- function(db, var.name, est.names, std.names, flag.na = TRUE) {
  stopifnot(all.equal(length(est.names), length(std.names)))
  score.names <- c("ME", "MAE", "RMSE", "NMSE", "LogS")
  score <- matrix(rep(-9999, length(score.names) * length(est.names)),
                  nrow = length(score.names), ncol = length(est.names))
  rownames(score) <- score.names
  colnames(score) <- est.names
  
  Z <- db.extract(db, names = var.name, flag.compress = TRUE)
  for (i in 1:length(est.names)) {
    Z.est <- db.extract(db, names = est.names[i], flag.compress = TRUE)
    Z.std <- db.extract(db, names = std.names[i], flag.compress = TRUE)
    score["ME", i]   <- mean(Z.est - Z, na.rm = flag.na)
    score["MAE", i]  <- mean(abs(Z.est - Z), na.rm = flag.na)
    score["RMSE", i] <- sqrt(mean((Z.est - Z) ^ 2, na.rm = flag.na))
    score["NMSE", i] <- mean(((Z.est - Z) / Z.std) ^ 2, na.rm = flag.na)
    score["LogS", i] <- mean(1/2 * log(2 * pi * Z.est ^ 2) + 1/2 * ((Z.est - Z) / Z.std) ^ 2, na.rm = flag.na)
  }
  score
}


## ----data-nocorr-homo, echo=FALSE, output="hide"------------------------------
kable(rbind(
  db.stat.multi(db.locate(data, names = "sel_data", loctype = "sel"),
                names = c("V1", "V2", "Z")),
  db.stat.multi(db.locate(data, names = "sel_target", loctype = "sel"),
                names = c("V1", "V2", "Z"))), 
  caption = "Statistiques des données initiales sans correlation et homotopiques")


## ----estim-nocorr-homo, , echo=TRUE, ouput='hide'-----------------------------
est.k <- ""
data <- db.delete(data, names = c("ok*", "ck*"))
data <- QC_LC_V12(cas = est.k)

# Score des estimations
Z.nm <- paste0("Z", est.k)
kable(db.score(data, 
               var.name = "Z", 
               est.names = paste(c("ck", "ok"), Z.nm, "estim",sep ="."),
               std.names = paste(c("ck", "ok"), Z.nm, "stdev",sep ="."),
               flag.na = TRUE), 
      caption = "Scores d'estimation pour le cas homotopique sans corrélation")


## ----data-nocorr-hetero, echo=FALSE, output="hide"----------------------------
kable(rbind(
  db.stat.multi(db.locate(data, names = "sel_data", loctype = "sel"),
                names = paste0(c("V1", "V2", "Z"), "_obs")),
  db.stat.multi(db.locate(data, names = "sel_target", loctype = "sel"),
                names = c("V1", "V2", "Z"))), 
  caption = "Statistiques des données initiales sans correlation et heterotopiques")


## ----estim-nocorr-hetero, echo=FALSE, ouput='hide'----------------------------
est.k <- "_obs"
data <- db.delete(data, names = c("ok*", "ck*"))
data <- QC_LC_V12(cas = est.k)
# Score des estimations
Z.nm <- paste0("Z", est.k)
kable(db.score(data, 
               var.name = "Z", 
               est.names = paste(c("ck", "ok"), Z.nm, "estim", sep="."),
               std.names = paste(c("ck", "ok"), Z.nm, "stdev", sep="."),
               flag.na = TRUE), 
      caption = "Scores d'estimation pour le cas hétérotopique sans corrélation")


## ----creating-data-corr, echo=FALSE, results='hide'---------------------------
data <- db.delete(data, names = c("Z", "Z_obs", "ok*", "ck*", "V*_obs"))
data <- db.locerase(data, loctype = "sel")
# corrélation entre les composantes de bases
data.rho <- 0.7
data.lu <- matrix(c(1, data.rho, 0, sqrt(1 - data.rho ^ 2)), 
                  nrow = 2, ncol = 2)
rownames(data.lu) <- c("V1", "V2")
colnames(data.lu) <- c("V1", "V2")
model.V <- model.CL(model = model.V, mat.CL = data.lu)
# données homotopiques
data <- db.add.CL(data, Y.names = colnames(data.lu), 
                  Z.names = paste0("VV", 1:2), mat.CL = data.lu)
data <- db.delete(data, names = paste0("V", 1:2))
data <- db.rename(data, name = "VV1", newname = "V1")
data <- db.rename(data, name = "VV2", newname = "V2")
data <- db.add.CL(data, Y.names = colnames(data.cl), 
                  Z.names = rownames(data.cl), mat.CL = data.cl)
# données hétérotopiques
for (i in 1:2) {
  W_obs <- db.extract(data, names = paste0("V", i), 
                      flag.compress = FALSE)
  W_obs[sample(x = c(TRUE, FALSE), size = data$nech,  replace = TRUE,
               prob = c(1/2, 1/2))] <- NA
  data <- db.add(data, W_obs)
  data <- db.rename(data, name = "W_obs", 
                    newname = paste0("V", i, "_obs"))
  }
data <- db.add.CL(data, Y.names = paste0("V", 1:2, "_obs"), 
                  Z.names = "Z_obs", mat.CL = data.cl)
data <- db.locerase(data, loctype = "z")
rm(W_obs)


## ----data-corr-homo, echo=FALSE, output='hide'--------------------------------
kable(rbind(
  db.stat.multi(db.locate(data, names = "sel_data", loctype = "sel"),
                names = c("V1", "V2", "Z")),
  db.stat.multi(db.locate(data, names = "sel_target", loctype = "sel"),
                names = c("V1", "V2", "Z")) ), 
  caption = "Statistiques des données initiales avec correlation et homotopiques")


## ----estim-corr-homo, echo=FALSE, ouput='hide'--------------------------------
est.k <- "_obs"
data <- db.delete(data, names = c("ok*", "ck*"))
data <- QC_LC_V12(cas = est.k)
# Score des estimations
Z.nm <- paste0("Z", est.k)
kable(db.score(data, 
               var.name = "Z", 
               est.names = paste(c("ck", "ok"), Z.nm, "estim", sep="."),
               std.names = paste(c("ck", "ok"), Z.nm, "stdev", sep="."),
               flag.na = TRUE), 
      caption = "Scores d'estimation pour le cas homotopique avec corrélation")


## ----data-corr-hetero, echo=FALSE, output='hide'------------------------------
kable(rbind(
  db.stat.multi(db.locate(data, names = "sel_data", loctype = "sel"),
                names = paste0(c("V1", "V2", "Z"), "_obs")),
  db.stat.multi(db.locate(data, names = "sel_target", loctype = "sel"),
                names = c("V1", "V2", "Z"))), 
  caption = "Statistiques des données initiales avec correlation et hétérotopiques")


## ----estim-corr-hetero, echo=FALSE, output="hide"-----------------------------
est.k <- "_obs"
data <- db.delete(data, names = c("ok*", "ck*"))
data <- QC_LC_V12(cas = est.k)
# Score des estimations
Z.nm <- paste0("Z", est.k)
kable(db.score(data, 
               var.name = "Z", 
               est.names = paste(c("ck", "ok"), Z.nm, "estim", sep="."),
               std.names = paste(c("ck", "ok"), Z.nm, "stdev", sep="."),
               flag.na = TRUE), 
      caption = "Scores d'estimation pour le cas hétérotopique avec corrélation")