paket <- c("sp", "sf", "gstat", "terra", "ggplot2", "viridis", "patchwork", "knitr")
baru  <- paket[!paket %in% rownames(installed.packages())]
if (length(baru) > 0) install.packages(baru)

library(sp)
library(sf)
library(gstat)
library(terra)
library(ggplot2)
library(viridis)
library(patchwork)
library(knitr)

# 1. Closest observation: nilai titik observasi terdekat.
pred_closest <- function(train, newdata) {
  idx <- st_nearest_feature(newdata, train)
  train$vble[idx]
}

# 2. IDW: seluruh titik dipakai, bobot = 1 / jarak^p (p = 2 sesuai materi)
pred_idw <- function(train, newdata, p = 2) {
  g <- gstat(formula = vble ~ 1, locations = train,
             nmax = nrow(train), set = list(idp = p))
  predict(g, newdata, debug.level = 0)$var1.pred
}

# 3. Nearest neighbors: rata-rata k tetangga terdekat (idp = 0 -> bobot sama)
pred_nn <- function(train, newdata, k = 5) {
  g <- gstat(formula = vble ~ 1, locations = train,
             nmax = k, set = list(idp = 0))
  predict(g, newdata, debug.level = 0)$var1.pred
}

# 4. Ordinary kriging: variogram dihitung dan di-fit ulang dari data train.
fit_vgm <- function(train) {
  bb    <- st_bbox(train)
  diagn <- sqrt((bb[3] - bb[1])^2 + (bb[4] - bb[2])^2)   # diagonal wilayah (m)
  vv    <- var(train$vble)
  vinit <- vgm(psill = 0.7 * vv, model = "Sph",
               range = as.numeric(diagn) / 4, nugget = 0.3 * vv)
  v0 <- variogram(vble ~ 1, train)
  fv <- tryCatch(suppressWarnings(fit.variogram(v0, model = vinit)),
                 error = function(e) vinit)
  # cadangan jika fitting gagal atau menghasilkan parameter tidak valid
  if (isTRUE(attr(fv, "singular")) || any(fv$psill < 0) || fv$range[2] <= 0) fv <- vinit
  fv
}

pred_ok <- function(train, newdata) {
  fv <- fit_vgm(train)
  krige(vble ~ 1, locations = train, newdata = newdata,
        model = fv, debug.level = 0)$var1.pred
}

RMSE <- function(obs, pred) sqrt(mean((obs - pred)^2))

K  <- 5
n  <- nrow(d)
kf <- sample(rep(1:K, length.out = n))   # penetapan lipatan acak

hasil <- data.frame(fold = 1:K, Closest = NA, NN = NA, IDW = NA,
                    Kriging = NA, Ens_Sama = NA, Ens_InvRMSE = NA)

for (k in 1:K) {
  test  <- d[kf == k, ]
  train <- d[kf != k, ]

  p1 <- pred_closest(train, test)
  p2 <- pred_nn(train, test)
  p3 <- pred_idw(train, test)
  p4 <- pred_ok(train, test)

  r <- c(RMSE(test$vble, p1), RMSE(test$vble, p2),
         RMSE(test$vble, p3), RMSE(test$vble, p4))

  # Ensemble bobot sama
  p5 <- (p1 + p2 + p3 + p4) / 4
  # Ensemble bobot 1/RMSE
  w  <- (1 / r) / sum(1 / r)
  p6 <- w[1] * p1 + w[2] * p2 + w[3] * p3 + w[4] * p4

  hasil[k, -1] <- c(r, RMSE(test$vble, p5), RMSE(test$vble, p6))
}

kable(hasil, digits = 4, caption = "RMSE per lipatan (skala log curah hujan harian)")

ringkas <- data.frame(
  Kode   = c("Closest", "NN", "IDW", "Kriging", "Ens_Sama", "Ens_InvRMSE"),
  Metode = c("Closest Observation", "Nearest Neighbors (k=5)", "IDW (p=2)",
             "Ordinary Kriging", "Ensemble (bobot sama)", "Ensemble (bobot 1/RMSE)"),
  RMSE_rata2 = colMeans(hasil[, -1]),
  stringsAsFactors = FALSE
)
ringkas <- ringkas[order(ringkas$RMSE_rata2), ]
rownames(ringkas) <- NULL
kable(ringkas[, c("Metode", "RMSE_rata2")], digits = 4,
      caption = "Rata-rata RMSE 5-fold CV (diurutkan, terkecil = terbaik)")

rg <- ringkas
rg$Metode <- factor(rg$Metode, levels = rev(rg$Metode))
ggplot(rg, aes(x = Metode, y = RMSE_rata2, fill = RMSE_rata2)) +
  geom_col(width = 0.7) +
  geom_text(aes(label = round(RMSE_rata2, 3)), hjust = -0.1) +
  coord_flip(ylim = c(0, max(rg$RMSE_rata2) * 1.15)) +
  scale_fill_viridis_c(direction = -1, guide = "none") +
  labs(title = "Perbandingan RMSE (5-fold CV)", x = NULL,
       y = "RMSE (log curah hujan harian)") +
  theme_bw()

# Bobot ensemble untuk peta akhir: 1/RMSE rata-rata hasil CV (tanpa mengintip data uji)
rm_cv <- colMeans(hasil[, c("Closest", "NN", "IDW", "Kriging")])
w_cv  <- (1 / rm_cv) / sum(1 / rm_cv)
round(w_cv, 3)

p_ens <- w_cv[["Closest"]] * p_closest + w_cv[["NN"]] * p_nn +
         w_cv[["IDW"]] * p_idw + w_cv[["Kriging"]] * p_ok

lim <- range(c(p_closest, p_nn, p_idw, p_ok, p_ens))
xy  <- st_coordinates(coop)

peta <- function(pred, judul) {
  df <- data.frame(X = xy[, 1], Y = xy[, 2], pred = as.numeric(pred))
  ggplot(df, aes(X, Y, fill = pred)) +
    geom_raster() +
    scale_fill_viridis_c(name = "log(CH harian)", limits = lim) +
    coord_equal() +
    labs(title = judul, x = NULL, y = NULL) +
    theme_bw(base_size = 9)
}

(peta(p_closest, "Closest Observation") | peta(p_nn, "Nearest Neighbors (k=5)")) /
(peta(p_idw, "IDW (p = 2)")             | peta(p_ok, "Ordinary Kriging")) /
(peta(p_ens, "Ensemble (bobot 1/RMSE CV)") | plot_spacer()) +
  plot_layout(guides = "collect")