library(sp)
library(sf)
library(terra)
library(gstat)
library(dplyr)
library(ggplot2)
library(viridis)

data(meuse)
data(meuse.grid)

# titik pengamatan, vble = log(zinc)
d <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992) %>%
  transmute(vble = log(zinc))

# lokasi yang akan diprediksi (grid 40 m x 40 m)
coop <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992, remove = FALSE)

# batas wilayah untuk diagram Voronoi
map <- st_sf(geometry = st_buffer(
  st_convex_hull(st_union(c(st_geometry(d), st_geometry(coop)))), 100))

c(jumlah_titik = nrow(d), jumlah_lokasi_prediksi = nrow(coop))

summary(d$vble)

ggplot() +
  geom_tile(data = meuse.grid, aes(x, y), fill = "grey90", width = 40, height = 40) +
  geom_sf(data = d, aes(color = vble), size = 2) +
  scale_color_viridis(name = "log(zinc)") +
  theme_bw() + labs(title = "Sebaran titik pengamatan log(zinc)")

RMSE <- function(observed, predicted) sqrt(mean((observed - predicted)^2))

# 1. Closest Observation: ambil nilai titik terdekat (poligon Voronoi)
closest_obs <- function(train, test) {
  v <- st_as_sf(terra::voronoi(terra::vect(train), bnd = terra::vect(map)))
  st_crs(v) <- st_crs(train)
  p <- st_join(test, v["vble"], join = st_within, suffix = c("", ".v"))
  p[[if ("vble.v" %in% names(p)) "vble.v" else "vble"]]
}

# 2. IDW: semua titik ikut, pangkat jarak (idp) = 1
idw_p <- function(train, test) {
  g <- gstat(formula = vble ~ 1, locations = train, nmax = nrow(train),
             set = list(idp = 1, debug = 0))
  predict(g, test)$var1.pred
}

# 3. Nearest Neighbors: rata-rata 5 titik terdekat, bobot sama (idp = 0)
nn_p <- function(train, test) {
  g <- gstat(formula = vble ~ 1, locations = train, nmax = 5,
             set = list(idp = 0, debug = 0))
  predict(g, test)$var1.pred
}

# 4. Ordinary Kriging: variogram spherical, nilai awal sama seperti di modul
fit_vgm <- function(train) {
  v <- variogram(vble ~ 1, train)
  fit.variogram(v, vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1))
}
kr_p <- function(train, test) {
  g <- gstat(formula = vble ~ 1, data = train, model = fit_vgm(train),
             set = list(debug = 0))
  predict(g, test)$var1.pred
}

v <- st_as_sf(terra::voronoi(terra::vect(d), bnd = terra::vect(map)))
st_crs(v) <- st_crs(d)

ggplot() +
  geom_sf(data = v, aes(fill = vble), alpha = 0.7, color = "white", linewidth = 0.1) +
  geom_sf(data = d, size = 0.6) +
  scale_fill_viridis(name = "log(zinc)") +
  theme_bw() + labs(title = "Diagram Voronoi")

metode <- function(train, test) {
  data.frame(closest = closest_obs(train, test), idw = idw_p(train, test),
             nn = nn_p(train, test), krig = kr_p(train, test))
}

# bobot invers RMSE dari 5-fold CV di dalam data latih
bobot_cv <- function(train, k = 5) {
  kf <- dismo::kfold(nrow(train), k)
  hasil <- do.call(rbind, lapply(1:k, function(i) {
    te <- train[kf == i, ]
    cbind(obs = te$vble, metode(train[kf != i, ], te))
  }))
  r <- sapply(hasil[, -1], function(p) RMSE(hasil$obs, p))
  w <- 1 / r
  w / sum(w)
}

set.seed(123)
w_all <- bobot_cv(d)

P <- metode(d, coop)
P$ens_sama  <- (P$closest + P$idw + P$nn) / 3
P$ens_bobot <- as.numeric(as.matrix(P[, c("closest", "idw", "nn", "krig")]) %*% w_all)

label <- c(closest = "Closest Observation", idw = "IDW", nn = "NN (k = 5)",
           ens_sama = "Ensembel A (bobot sama)", krig = "Kriging",
           ens_bobot = "Ensembel B (bobot RMSE)")
urut <- c("closest", "idw", "nn", "ens_sama", "krig", "ens_bobot")

peta <- do.call(rbind, lapply(urut, function(m)
  data.frame(x = coop$x, y = coop$y, pred = P[[m]], metode = label[m])))
peta$metode <- factor(peta$metode, levels = label[urut])

ggplot(peta, aes(x, y, fill = pred)) + geom_tile(width = 40, height = 40) +
  facet_wrap(~metode, ncol = 3) + coord_equal() +
  scale_fill_viridis(name = "log(zinc)") + theme_bw() +
  theme(axis.text = element_blank(), axis.ticks = element_blank(),
        axis.title = element_blank(), strip.text = element_text(size = 8))

knitr::kable(data.frame(Metode = c("Closest Observation", "IDW", "NN", "Kriging"),
                        Bobot = round(w_all, 3)), row.names = FALSE)

cv_run <- function(seed, k = 5) {
  set.seed(seed)
  kf <- dismo::kfold(nrow(d), k)
  do.call(rbind, lapply(1:k, function(i) {
    tr <- d[kf != i, ]; te <- d[kf == i, ]
    p <- metode(tr, te)
    w <- bobot_cv(tr)
    p$ens_sama  <- (p$closest + p$idw + p$nn) / 3
    p$ens_bobot <- as.numeric(as.matrix(p[, c("closest", "idw", "nn", "krig")]) %*% w)
    r <- sapply(p[urut], function(x) RMSE(te$vble, x))
    data.frame(fold = i, t(r))
  }))
}

hasil <- cv_run(123)
names(hasil)[-1] <- label[urut]
knitr::kable(hasil, digits = 4, row.names = FALSE, caption = "RMSE tiap lipatan (skala log)")

rata <- colMeans(hasil[, -1])
tab <- data.frame(Metode = names(rata), RMSE_rata_rata = round(rata, 4),
                  Peringkat = rank(rata), row.names = NULL)
knitr::kable(tab[order(tab$Peringkat), ], row.names = FALSE,
             caption = "Rata-rata RMSE 5 lipatan")

tab$Metode <- factor(tab$Metode, levels = tab$Metode[order(tab$RMSE_rata_rata)])
ggplot(tab, aes(Metode, RMSE_rata_rata, fill = Metode)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = sprintf("%.3f", RMSE_rata_rata)), vjust = -0.4) +
  scale_fill_viridis_d() + theme_bw() +
  labs(y = "RMSE rata-rata", x = NULL) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))