# COSC/MATH 201 - Empirical Models in Practice (Module 8.3)
# Choose a shape, fit it, then test it on data the model hasn't seen.
# Save boxbod.csv in the same folder as this script; pressure is built into R.

par(mfrow = c(2, 3), mar = c(2, 2, 2, 1))
curve(1 + 2 * x, 0, 5, main = "Linear", lwd = 2, col = "steelblue")
curve(x^2.5, 0, 5, main = "Power (p = 2.5)", lwd = 2, col = "steelblue")
curve(exp(0.8 * x), 0, 5, main = "Exponential", lwd = 2, col = "steelblue")
curve(10 * (1 - exp(-1.2 * x)), 0, 5, main = "Leveling off", lwd = 2,
      col = "steelblue")
curve(10 / (1 + 50 * exp(-2 * x)), 0, 5, main = "Logistic", lwd = 2,
      col = "steelblue")
par(mfrow = c(1, 1))

boxbod <- read.csv("boxbod.csv")
boxbod

rates <- seq(0.1, 1.5, by = 0.05)
ssr <- numeric(length(rates))

for (i in seq_along(rates)) {
  z <- 1 - exp(-rates[i] * boxbod$days)
  fit <- lm(boxbod$bod ~ 0 + z)
  ssr[i] <- sum(resid(fit)^2)
}

best_k <- rates[which.min(ssr)]
best_k

plot(rates, ssr, type = "o", pch = 16, col = "purple",
     xlab = "k (per day)", ylab = "Sum of squared residuals",
     main = "Searching for the best k")
abline(v = best_k, lty = 2)

boxbod$z <- 1 - exp(-best_k * boxbod$days)
level_fit <- lm(bod ~ 0 + z, data = boxbod)
best_a <- coef(level_fit)
best_a

day_grid <- seq(0, 20, by = 0.1)
plot(boxbod$days, boxbod$bod, pch = 16, cex = 1.3, col = "purple",
     xlim = c(0, 20), ylim = c(0, 250),
     xlab = "Incubation time (days)", ylab = "BOD (mg/L)",
     main = "Leveling-off model for BOD")
lines(day_grid, best_a * (1 - exp(-best_k * day_grid)), col = "darkgreen", lwd = 2)
abline(v = 10, lty = 3, col = "gray50")   # last measurement

new_days <- c(8, 12, 20)
data.frame(days = new_days, predicted_bod = round(best_a * (1 - exp(-best_k * new_days))))

nls(bod ~ a * (1 - exp(-k * days)), data = boxbod,
    start = list(a = 200, k = 0.5))

head(pressure)
plot(pressure$temperature, pressure$pressure, pch = 16, col = "steelblue",
     xlab = "Temperature (°C)", ylab = "Vapor pressure (mm Hg)",
     main = "Mercury vapor pressure")

pressure$kelvin <- pressure$temperature + 273.15
train <- pressure[pressure$temperature <= 260, ]
test <- pressure[pressure$temperature > 260, ]
nrow(train)
nrow(test)

expo_fit <- lm(log(pressure) ~ temperature, data = train)
cubic_fit <- lm(pressure ~ temperature + I(temperature^2) + I(temperature^3),
                data = train)
physics_fit <- lm(log(pressure) ~ I(1 / kelvin), data = train)

predictions <- data.frame(
  temperature = test$temperature,
  measured = test$pressure,
  exponential = round(exp(predict(expo_fit, test))),
  cubic = round(predict(cubic_fit, test)),
  physics = round(exp(predict(physics_fit, test)))
)
predictions

temp_grid <- seq(0, 360, by = 2)
grid <- data.frame(temperature = temp_grid, kelvin = temp_grid + 273.15)
plot(pressure$temperature, pressure$pressure, pch = 16,
     col = ifelse(pressure$temperature <= 260, "steelblue", "firebrick"),
     ylim = c(0, 1000), xlab = "Temperature (°C)",
     ylab = "Vapor pressure (mm Hg)", main = "Trained on blue points, tested on red")
lines(temp_grid, exp(predict(expo_fit, grid)), col = "orange", lwd = 2)
lines(temp_grid, predict(cubic_fit, grid), col = "gray40", lwd = 2, lty = 2)
lines(temp_grid, exp(predict(physics_fit, grid)), col = "darkgreen", lwd = 2)
abline(v = 270, lty = 3)
legend("topleft", c("Exponential", "Cubic", "Physics-based"),
       col = c("orange", "gray40", "darkgreen"), lwd = 2, lty = c(1, 2, 1),
       bty = "n")

growth_fraction <- 0.1
end_time <- 100
exact <- 100 * exp(growth_fraction * end_time)

simulate <- function(method, delta_t) {
  population <- 100
  for (step in 1:round(end_time / delta_t)) {
    if (method == "Euler") {
      population <- population + growth_fraction * population * delta_t
    } else if (method == "RK2") {
      change_1 <- growth_fraction * population * delta_t
      change_2 <- growth_fraction * (population + change_1) * delta_t
      population <- population + (change_1 + change_2) / 2
    } else {
      change_1 <- growth_fraction * population * delta_t
      change_2 <- growth_fraction * (population + change_1 / 2) * delta_t
      change_3 <- growth_fraction * (population + change_2 / 2) * delta_t
      change_4 <- growth_fraction * (population + change_3) * delta_t
      population <- population + (change_1 + 2 * change_2 + 2 * change_3 + change_4) / 6
    }
  }
  population
}

steps <- c(0.1, 0.2, 0.25, 0.5, 1)
methods <- c("Euler", "RK2", "RK4")
errors <- matrix(0, nrow = length(steps), ncol = length(methods),
                 dimnames = list(steps, methods))

for (m in methods) {
  for (i in seq_along(steps)) {
    errors[i, m] <- abs(simulate(m, steps[i]) - exact) / exact
  }
}
signif(errors, 3)

for (m in methods) {
  slope <- coef(lm(log(errors[, m]) ~ log(steps)))[2]
  cat(m, "error shrinks like delta_t to the power", round(slope, 2), "\n")
}

plot(log(steps), log(errors[, "Euler"]), type = "o", pch = 16, col = "firebrick",
     ylim = range(log(errors)), xlab = "ln(delta t)", ylab = "ln(relative error)",
     main = "Error vs. time step (log-log scale)")
lines(log(steps), log(errors[, "RK2"]), type = "o", pch = 16, col = "steelblue")
lines(log(steps), log(errors[, "RK4"]), type = "o", pch = 16, col = "darkgreen")
legend("bottomright", methods, col = c("firebrick", "steelblue", "darkgreen"),
       pch = 16, lwd = 1, bty = "n")