# Load packages
library(nlme); library(mgcv); library(nnet); library(cluster)

# Config
MAX_TORQUE_NM <- 450
impl <- c(F1="Plough", F2="Plough", F3="Plough", F4="Plough", F5="Plough",
          F6="Plough", F7="RotaryTiller", F8="Plough", F9="Plough", F10="RotaryTiller")

# Read & clean
read_f <- function(p) {
  d <- read.delim(p, sep="\t")
  d$Torque_Nm <- d$ActualEngine_PercTorque / 100 * MAX_TORQUE_NM
  d$HtcDraftSens1Perc[d$HtcDraftSens1Perc < -50 | d$HtcDraftSens1Perc > 150] <- NA
  d$HtcDraftSens2Perc[d$HtcDraftSens2Perc < -50 | d$HtcDraftSens2Perc > 150] <- NA
  d$HtcPositionSensPerc[d$HtcPositionSensPerc < -10 | d$HtcPositionSensPerc > 110] <- NA
  d$EngOilPress[d$EngOilPress <= 0] <- NA
  d$t_s <- d$time_ms / 1000
  d
}
files <- paste0("FIELD_", 1:10, "_TUMOSAN_81110P_TILLAGE_OP_RAW.txt")
dat <- lapply(file.path("./FIELD_FILES", files), read_f)
names(dat) <- paste0("F", 1:10)
for (nm in names(dat)) dat[[nm]]$field <- nm

# Hysteresis segmentation
seg_hyst <- function(d, min_dwell=3) {
  hp <- d$HtcPositionSensPerc
  hp[is.na(hp)] <- median(hp, na.rm=TRUE)
  lo <- quantile(hp, 0.25); hi <- quantile(hp, 0.60)
  st <- logical(length(hp))
  cur <- hp[1] <= lo
  for (i in seq_along(hp)) {
    if (hp[i] <= lo) cur <- TRUE
    if (hp[i] >= hi) cur <- FALSE
    st[i] <- cur
  }
  rl <- rle(st); v <- rl$values; l <- rl$lengths
  min_s <- min_dwell * 10
  repeat {
    ch <- FALSE
    for (i in seq_along(l)) {
      if (l[i] < min_s && length(l) > 1) {
        if (i == 1) { l[2] <- l[2]+l[1]; v <- v[-1]; l <- l[-1] }
        else if (i == length(l)) { l[i-1] <- l[i-1]+l[i]; v <- v[-i]; l <- l[-i] }
        else { l[i-1] <- l[i-1]+l[i]+l[i+1]; v <- v[-c(i,i+1)]; l <- l[-c(i,i+1)] }
        ch <- TRUE; break
      }
    }
    if (!ch || length(l) <= 1) break
  }
  rep(v, l)
}
for (nm in names(dat)) dat[[nm]]$is_work <- seg_hyst(dat[[nm]])

# Work-phase dataset
work <- do.call(rbind, lapply(names(dat), function(nm) {
  d <- dat[[nm]]
  data.frame(field=nm, impl=impl[nm], torque=d$ActualEngine_PercTorque[d$is_work],
             fuel=d$EngFuelRate[d$is_work], rpm=d$EngineSpeed[d$is_work],
             speed=d$WheelBasedVehicleSpeed[d$is_work],
             draft1=d$HtcDraftSens1Perc[d$is_work], draft2=d$HtcDraftSens2Perc[d$is_work])
}))
work <- work[complete.cases(work$torque),]
work$field <- factor(work$field)
work$power <- work$torque/100 * work$rpm

# ---- Mixed model (implement effect) ----
m <- lme(torque ~ impl, random=~1|field, data=work)
summary(m)$tTable
vc <- as.numeric(VarCorr(m)[,"Variance"])
icc <- vc[1]/sum(vc); cat("ICC:", round(icc,3), "\n")

# Bootstrap on field medians (Cliff's delta & CI)
meds <- aggregate(torque ~ field + impl, data=work, median)
cliff <- function(x,y) sum(outer(x, y, function(a,b) sign(a-b))) / (length(x)*length(y))
cat("Cliff's delta:", round(cliff(meds$torque[meds$impl=="Plough"], meds$torque[meds$impl=="RotaryTiller"]),3), "\n")

set.seed(123)
boot_diff <- replicate(2000, {
  s <- sample(levels(work$field), replace=TRUE)
  pm <- if(any(impl[s]=="Plough")) median(work$torque[work$field %in% s[impl[s]=="Plough"]]) else NA
  rm <- if(any(impl[s]=="RotaryTiller")) median(work$torque[work$field %in% s[impl[s]=="RotaryTiller"]]) else NA
  pm - rm
})
cat("Boot CI median diff:", round(quantile(boot_diff, c(.025,.975), na.rm=TRUE),1), "\n")

# ICC bootstrap
set.seed(99)
icc_boot <- replicate(500, {
  s <- sample(levels(work$field), replace=TRUE)
  tmp <- do.call(rbind, lapply(seq_along(s), function(i) data.frame(field=paste0("r",i), torque=work$torque[work$field==s[i]])))
  mm <- tryCatch(lme(torque~1, random=~1|field, data=tmp), error=function(e) NULL)
  if(!is.null(mm)) { v <- as.numeric(VarCorr(mm)[,"Variance"]); v[1]/sum(v) } else NA
})
cat("ICC 95% CI:", round(quantile(icc_boot, c(.025,.975), na.rm=TRUE),3), "\n")

# ---- GAM fuel model ----
m_lin <- lm(fuel ~ power, data=work)
m_gam <- bam(fuel ~ s(torque, k=8) + s(rpm, k=8) + s(speed, k=8) + impl, data=work, discrete=TRUE)
cat("Linear R2:", round(summary(m_lin)$r.squared,3), " GAM dev:", round(summary(m_gam)$dev.expl,3), "\n")

png("gam_terms.png", width=1400, height=500, res=150)
par(mfrow=c(1,3)); plot(m_gam, select=1, shade=TRUE); plot(m_gam, select=2, shade=TRUE); plot(m_gam, select=3, shade=TRUE)
dev.off()
# ---- AR(1) robustness check (Section 3.3/4.3, Table 2) ----
# VERIFIED: reproduces SE~10.4 baseline at both resolutions,
# AR1-corrected: 10s -> SE=2.73, p=0.0155, deltaAIC=-3797
#               5s  -> SE=2.15, p=0.0038, deltaAIC=-9104
agg_seq <- function(d, bin_s) {
  w <- d[d$is_work,]; w <- w[complete.cases(w$ActualEngine_PercTorque),]
  n_per_bin <- bin_s * 10
  grp <- ceiling(seq_len(nrow(w)) / n_per_bin)
  ag <- aggregate(ActualEngine_PercTorque ~ grp, data=w, mean)
  data.frame(field=d$field[1], impl=impl[d$field[1]], seqidx=ag$grp, torque=ag$ActualEngine_PercTorque)
}
run_ar1 <- function(bin_s) {
  agdat <- do.call(rbind, lapply(dat, agg_seq, bin_s=bin_s))
  agdat$field <- factor(agdat$field)
  agdat <- agdat[order(agdat$field, agdat$seqidx),]
  m0 <- lme(torque ~ impl, random=~1|field, data=agdat)
  m1 <- lme(torque ~ impl, random=~1|field, correlation=corAR1(form=~seqidx|field), data=agdat)
  cat("bin=", bin_s, "s  no-AR1 SE=", round(summary(m0)$tTable[2,2],2),
      " AR1 SE=", round(summary(m1)$tTable[2,2],2),
      " AR1 p=", round(summary(m1)$tTable[2,5],4),
      " deltaAIC=", round(AIC(m1)-AIC(m0)), "\n")
  list(m0=m0, m1=m1)
}
r10 <- run_ar1(10)
r5  <- run_ar1(5)

# ---- PCA + k-means ----
work$draft <- rowMeans(cbind(work$draft1, work$draft2), na.rm=TRUE)
pca_vars <- c("torque","rpm","fuel","speed","draft")
pca_df <- work[complete.cases(work[,pca_vars]), pca_vars]
pca <- prcomp(scale(pca_df))
set.seed(1)
km <- kmeans(pca$x[,1:2], centers=2, nstart=25)
work$cluster <- NA
work[complete.cases(work[,pca_vars]), "cluster"] <- km$cluster
prop.table(table(cluster=work$cluster, impl=work$impl), 1)

png("pca_clusters.png", width=1600, height=800, res=150)
idx <- sample(nrow(pca$x), 20000)
plot(pca$x[idx,1], pca$x[idx,2], col=adjustcolor(km$cluster[idx]+1, 0.3), pch=16, cex=0.4)
dev.off()

# ---- LOOCV neural network (lagged forecasting) ----
agg <- function(d) {
  grp <- ceiling((1:nrow(d))/10)
  data.frame(torque=tapply(d$ActualEngine_PercTorque, grp, mean),
             rpm=tapply(d$EngineSpeed, grp, mean),
             speed=tapply(d$WheelBasedVehicleSpeed, grp, mean),
             field=d$field[1])
}
dat_agg <- lapply(dat, agg)

make_ts <- function(d, k=5, h=3) {
  n <- nrow(d); idx <- k:(n-h)
  if(length(idx)<10) return(NULL)
  idx <- if(length(idx)>15000) sample(idx, 15000) else idx
  X <- t(sapply(idx, function(i) d$torque[i:(i-k+1)]))
  data.frame(X, rpm=d$rpm[idx], speed=d$speed[idx], y=d$torque[idx+h], field=d$field[1])
}
ts_list <- lapply(dat_agg, make_ts); ts_list <- ts_list[!sapply(ts_list, is.null)]

loocv <- lapply(names(ts_list), function(te) {
  tr <- do.call(rbind, ts_list[names(ts_list)!=te])
  te_df <- ts_list[[te]]
  mu <- colMeans(tr[,1:7]); sdv <- apply(tr[,1:7], 2, sd)
  norm <- function(x) sweep(sweep(x[,1:7], 2, mu, "-"), 2, sdv, "/")
  Xtr <- norm(tr); ytr <- tr$y
  Xte <- norm(te_df); yte <- te_df$y
  set.seed(42)
  nn <- nnet(Xtr, ytr/120, size=6, decay=0.005, linout=TRUE, maxit=200, trace=FALSE)
  pred_nn <- predict(nn, Xte)*120
  lm_vars <- setdiff(names(tr), c("y","field"))
pred_lm <- predict(lm(as.formula(paste("y ~", paste(lm_vars, collapse=" + "))), data=tr), te_df)
  pred_naive <- te_df[,1]
  data.frame(test=te, r2_naive=1-sum((pred_naive-yte)^2)/sum((yte-mean(yte))^2),
             r2_lm=1-sum((pred_lm-yte)^2)/sum((yte-mean(yte))^2),
             r2_nn=1-sum((pred_nn-yte)^2)/sum((yte-mean(yte))^2))
})
loocv <- do.call(rbind, loocv)
print(loocv)
cat("Mean NN R2:", round(mean(loocv$r2_nn),3), " SD:", round(sd(loocv$r2_nn),3), "\n")

# ---- Boxplots & duty cycles ----
png("torque_boxplot.png", width=1200, height=700, res=150)
boxplot(torque~field, data=work, col="orange", outline=FALSE)
dev.off()

png("duty_cycle.png", width=1400, height=900, res=150)
par(mfrow=c(5,2), mar=c(4,4,2,1))
for(nm in paste0("F",1:10)) {
  d <- dat[[nm]]; tq <- d$ActualEngine_PercTorque[d$is_work]
  tq <- tq[!is.na(tq) & tq>=0 & tq<=120]
  tab <- prop.table(table(cut(tq, seq(0,120,10), include.lowest=TRUE)))*100
  barplot(tab, main=paste(nm, impl[nm]), ylab="%")
}
dev.off()