write.csv(as.data.frame(rsm::decode.data(rsm::heli)), "heli.csv", row.names = FALSE)

library(rsm)

# Data, response, predictors and block
data_raw   <- as.data.frame(decode.data(heli))
response   <- "ave"
predictors <- c("A", "R", "W", "L")
block      <- "block"

# Find the CCD centre and factorial half-width
A <- sort(unique(na.omit(data_raw$A)), decreasing = TRUE)
R <- sort(unique(na.omit(data_raw$R)), decreasing = TRUE)
W <- sort(unique(na.omit(data_raw$W)), decreasing = TRUE)
L <- sort(unique(na.omit(data_raw$L)), decreasing = TRUE)

A_centre <- A[3]; A_half <- A[2] - A[3]
R_centre <- R[3]; R_half <- R[2] - R[3]
W_centre <- W[3]; W_half <- W[2] - W[3]
L_centre <- L[3]; L_half <- L[2] - L[3]

coded_data <- coded.data(
  data_raw,
  x1 ~ (A - A_centre) / A_half,
  x2 ~ (R - R_centre) / R_half,
  x3 ~ (W - W_centre) / W_half,
  x4 ~ (L - L_centre) / L_half
)

data_raw$block <- as.numeric(data_raw$block)   # to match the app's intercept

model_fit <- rsm(ave ~ block + SO(x1, x2, x3, x4), data = coded_data)
summary(model_fit)

# ---- Default view: partitioned ANOVA with lack of fit and pure error ----
summary(model_fit)$lof

# ---- Alternative view: term-wise sequential (Type I) ANOVA ----
lm_fit <- lm(ave ~ block + x1 + x2 + x3 + x4 +
               x1:x2 + x1:x3 + x1:x4 + x2:x3 + x2:x4 + x3:x4 +
               I(x1^2) + I(x2^2) + I(x3^2) + I(x4^2),
             data = coded_data)

anova_table <- anova(lm_fit)
anova_table

# dist = seq(-5, 5, by = 0.5) 
canonical_result <- rsm::canonical.path(model_fit)

library(rsm)
library(desirability)

df <- as.data.frame(decode.data(heli))

# Fitting a second order model
fit_ave   <- rsm(ave   ~ SO(A, R, W, L), data = df)
fit_logSD <- rsm(logSD ~ SO(A, R, W, L), data = df)

p_ave   <- predict(fit_ave,   df)
p_logSD <- predict(fit_logSD, df)

# Turn each response into a 0-1 desirability score.
# low and high values are adjusted to the default values
d_ave   <- dMax(low = min(p_ave)   - 0.1 * diff(range(p_ave)),
                high = max(p_ave)   + 0.1 * diff(range(p_ave)),   tol = 0.01)
d_logSD <- dMin(low = min(p_logSD) - 0.1 * diff(range(p_logSD)),
                high = max(p_logSD) + 0.1 * diff(range(p_logSD)), tol = 0.01)

# 3. Overall desirability D at any setting = geometric mean of the two scores
overall_D <- function(x) {
  new_row <- data.frame(A = x[1], R = x[2], W = x[3], L = x[4])
  sqrt(predict(d_ave,   predict(fit_ave,   new_row)) *
       predict(d_logSD, predict(fit_logSD, new_row)))
}

# 4. Search for the settings with the highest D, staying inside the tested range
lower_b <- c(min(df$A), min(df$R), min(df$W), min(df$L))
upper_b <- c(max(df$A), max(df$R), max(df$W), max(df$L))
start   <- c(mean(df$A), mean(df$R), mean(df$W), mean(df$L))

best <- optim(start, overall_D, method = "L-BFGS-B",
              lower = lower_b, upper = upper_b,
              control = list(fnscale = -1))   # fnscale = -1 means MAXIMISE

# 5. Results
xopt    <- best$par
new_row <- data.frame(A = xopt[1], R = xopt[2], W = xopt[3], L = xopt[4])

xopt                                      # optimum settings
predict(fit_ave,   new_row)               # predicted ave there
predict(fit_logSD, new_row)               # predicted logSD there
best$value