############################################################
# PRACTICAL FOR CLASSES 1–3
# Linear Programming in R (BASE R / R GUI / R console)
#
# Recommended R version:
# R 4.5.3
#
# Official download pages:
# Main CRAN page:
# https://cran.r-project.org/
#
# Windows:
# https://cran.r-project.org/bin/windows/base/
#
# macOS:
# https://cran.r-project.org/bin/macosx/
#
# How to run in BASE R (not RStudio):
# 1. Save this file, for example as lp_practical_base_R.R
# 2. Open R
# 3. Set your working directory if needed:
#    setwd("your_folder_path")
# 4. Run:
#    source("lp_practical_base_R.R")
#
# This script is written for standard CRAN R 4.5.3 and uses only:
# - ggplot2
# - lpSolve
#
# What this script covers
# 1. Modelling a decision problem
# 2. Feasible region and graphical solution
# 3. Binding constraints, slack, and shadow prices
# 4. A second example with non-zero slack
############################################################

############################
# 0. PACKAGE SETUP
############################

required_packages <- c("ggplot2", "lpSolve")

for (pkg in required_packages) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg, repos = "https://cloud.r-project.org")
    library(pkg, character.only = TRUE)
  }
}

############################
# 1. CLASS 1: MODELLING
############################

cat("\n==============================\n")
cat("CLASS 1: MODELLING A DECISION PROBLEM\n")
cat("==============================\n")

# Decision variables
# x1 = units of Product A
# x2 = units of Product B
#
# Objective
# Maximise profit:
# z = 40x1 + 30x2
#
# Constraints
# 2x1 + x2 <= 100   (machine time)
# x1 + 2x2 <= 80    (labour time)
# x1, x2 >= 0

objective <- c(40, 30)

constraint_matrix <- matrix(c(
  2, 1,
  1, 2
), nrow = 2, byrow = TRUE)

constraint_directions <- c("<=", "<=")
constraint_rhs <- c(100, 80)

cat("\nDecision variables:\n")
cat("x1 = units of Product A\n")
cat("x2 = units of Product B\n")

cat("\nObjective function:\n")
cat("Maximise z = 40x1 + 30x2\n")

cat("\nConstraints:\n")
cat("2x1 + x2 <= 100   (machine time)\n")
cat("x1 + 2x2 <= 80    (labour time)\n")
cat("x1 >= 0, x2 >= 0\n")

############################
# 2. CLASS 2: FEASIBLE REGION
############################

cat("\n==============================\n")
cat("CLASS 2: FEASIBLE REGION AND GRAPHICAL LP\n")
cat("==============================\n")

# A point is feasible if it satisfies every constraint.
is_feasible <- function(x1, x2) {
  (2 * x1 + x2 <= 100 + 1e-9) &&
    (x1 + 2 * x2 <= 80 + 1e-9) &&
    (x1 >= -1e-9) &&
    (x2 >= -1e-9)
}

# Generate a grid of points to visualise the feasible region.
grid <- expand.grid(
  x1 = seq(0, 60, by = 0.25),
  x2 = seq(0, 50, by = 0.25)
)

grid$feasible <- mapply(is_feasible, grid$x1, grid$x2)

cat("\nNumber of grid points tested:", nrow(grid), "\n")
cat("Number of feasible grid points:", sum(grid$feasible), "\n")

# Corner points:
# (0,0)
# (50,0) from 2x1 + x2 = 100 and x2 = 0
# (0,40) from x1 + 2x2 = 80 and x1 = 0
# (40,20) from solving:
# 2x1 + x2 = 100
# x1 + 2x2 = 80

corner_points <- data.frame(
  x1 = c(0, 50, 0, 40),
  x2 = c(0, 0, 40, 20)
)

corner_points$z <- 40 * corner_points$x1 + 30 * corner_points$x2

cat("\nCorner points and objective values:\n")
print(corner_points)

best_corner <- corner_points[which.max(corner_points$z), ]
cat("\nBest corner from graphical method:\n")
print(best_corner)

# Constraint lines for plotting
x <- seq(0, 60, length.out = 400)
line1 <- 100 - 2 * x
line2 <- (80 - x) / 2

# Plot 1: Feasible region
p1 <- ggplot() +
  geom_point(
    data = subset(grid, feasible),
    aes(x = x1, y = x2),
    alpha = 0.12,
    size = 0.5
  ) +
  geom_line(
    data = data.frame(x = x, y = line1),
    aes(x = x, y = y),
    linewidth = 1
  ) +
  geom_line(
    data = data.frame(x = x, y = line2),
    aes(x = x, y = y),
    linewidth = 1
  ) +
  geom_path(
    data = rbind(corner_points, corner_points[1, ]),
    aes(x = x1, y = x2),
    linewidth = 1
  ) +
  geom_point(
    data = corner_points,
    aes(x = x1, y = x2),
    size = 3
  ) +
  geom_text(
    data = corner_points,
    aes(
      x = x1,
      y = x2,
      label = paste0("(", x1, ", ", x2, ")")
    ),
    vjust = -1,
    size = 4
  ) +
  geom_point(
    data = best_corner,
    aes(x = x1, y = x2),
    size = 4
  ) +
  geom_text(
    data = best_corner,
    aes(
      x = x1,
      y = x2,
      label = paste0("Optimal: (", x1, ", ", x2, ")")
    ),
    nudge_x = 1,
    nudge_y = -3,
    size = 4
  ) +
  labs(
    title = "Feasible Region and Corner Points",
    x = "x1 (Product A)",
    y = "x2 (Product B)"
  ) +
  coord_cartesian(xlim = c(0, 60), ylim = c(0, 50)) +
  theme_minimal(base_size = 13)

print(p1)

# Plot 2: Objective-function lines
objective_levels <- c(1000, 1600, 2000, 2200)

objective_data <- do.call(
  rbind,
  lapply(objective_levels, function(k) {
    data.frame(
      x = x,
      y = (k - 40 * x) / 30,
      label = paste0("40x1 + 30x2 = ", k)
    )
  })
)

p2 <- ggplot() +
  geom_point(
    data = subset(grid, feasible),
    aes(x = x1, y = x2),
    alpha = 0.12,
    size = 0.5
  ) +
  geom_line(
    data = data.frame(x = x, y = line1),
    aes(x = x, y = y),
    linewidth = 1
  ) +
  geom_line(
    data = data.frame(x = x, y = line2),
    aes(x = x, y = y),
    linewidth = 1
  ) +
  geom_line(
    data = objective_data,
    aes(x = x, y = y, linetype = label),
    linewidth = 1
  ) +
  geom_point(
    data = best_corner,
    aes(x = x1, y = x2),
    size = 4
  ) +
  labs(
    title = "Objective Function Lines",
    subtitle = "Equal-profit lines move outward until they touch the feasible region",
    x = "x1 (Product A)",
    y = "x2 (Product B)",
    linetype = "Objective line"
  ) +
  coord_cartesian(xlim = c(0, 60), ylim = c(0, 50)) +
  theme_minimal(base_size = 13)

print(p2)

############################
# 3. CLASS 3: SOLVE AND INTERPRET
############################

cat("\n==============================\n")
cat("CLASS 3: SOLUTION, SLACK, SHADOW PRICES\n")
cat("==============================\n")

# Solve the LP using lpSolve
lp_result <- lpSolve::lp(
  direction = "max",
  objective.in = objective,
  const.mat = constraint_matrix,
  const.dir = constraint_directions,
  const.rhs = constraint_rhs,
  compute.sens = TRUE
)

cat("\nSolver status (0 means optimal):\n")
print(lp_result$status)

cat("\nOptimal objective value:\n")
print(lp_result$objval)

cat("\nOptimal decision variables:\n")
print(lp_result$solution)

used_resources <- as.numeric(constraint_matrix %*% lp_result$solution)
slack <- constraint_rhs - used_resources

results_table <- data.frame(
  Constraint = c("Machine time", "Labour time"),
  Capacity = constraint_rhs,
  Used = used_resources,
  Slack = slack
)

results_table$Binding <- ifelse(abs(results_table$Slack) < 1e-8, "Yes", "No")

cat("\nConstraint usage and slack:\n")
print(results_table)

dual_prices <- lp_result$duals[1:length(constraint_rhs)]

dual_table <- data.frame(
  Constraint = c("Machine time", "Labour time"),
  ShadowPrice = dual_prices
)

cat("\nShadow prices:\n")
print(dual_table)

cat("\nInterpretation of shadow prices:\n")
for (i in 1:nrow(dual_table)) {
  cat(
    paste0(
      "- ", dual_table$Constraint[i], ": if this resource increases by 1 unit, ",
      "the optimal profit changes by about ",
      round(dual_table$ShadowPrice[i], 4), ".\n"
    )
  )
}

# Plot 3: Resource usage vs capacity
resource_plot_data <- rbind(
  data.frame(Constraint = results_table$Constraint, Type = "Capacity", Value = results_table$Capacity),
  data.frame(Constraint = results_table$Constraint, Type = "Used", Value = results_table$Used)
)

p3 <- ggplot(resource_plot_data, aes(x = Constraint, y = Value, fill = Type)) +
  geom_col(position = "dodge") +
  geom_text(
    aes(label = round(Value, 1)),
    position = position_dodge(width = 0.9),
    vjust = -0.4,
    size = 4
  ) +
  labs(
    title = "Resource Usage vs Capacity",
    subtitle = "Both constraints are binding in this example",
    x = "Constraint",
    y = "Units"
  ) +
  theme_minimal(base_size = 13)

print(p3)

# Plot 4: Slack
# In this first example both slack values are zero.
p4 <- ggplot(results_table, aes(x = Constraint, y = Slack)) +
  geom_col() +
  geom_text(aes(label = round(Slack, 2)), vjust = -0.4, size = 4) +
  labs(
    title = "Slack by Constraint",
    subtitle = "Slack = unused amount of a resource (here both are zero)",
    x = "Constraint",
    y = "Slack"
  ) +
  theme_minimal(base_size = 13)

print(p4)

# Plot 5: Shadow prices
p5 <- ggplot(dual_table, aes(x = Constraint, y = ShadowPrice)) +
  geom_col() +
  geom_text(aes(label = round(ShadowPrice, 2)), vjust = -0.4, size = 4) +
  labs(
    title = "Shadow Prices",
    subtitle = "Value of one additional unit of each resource",
    x = "Constraint",
    y = "Shadow price"
  ) +
  theme_minimal(base_size = 13)

print(p5)

############################
# 4. SENSITIVITY: BEFORE / AFTER
############################

cat("\n==============================\n")
cat("SENSITIVITY: ORIGINAL VS RELAXED CONSTRAINT\n")
cat("==============================\n")

# Relax the first constraint:
# 2x1 + x2 <= 120 instead of 100
is_feasible_modified <- function(x1, x2) {
  (2 * x1 + x2 <= 120 + 1e-9) &&
    (x1 + 2 * x2 <= 80 + 1e-9) &&
    (x1 >= -1e-9) &&
    (x2 >= -1e-9)
}

grid_compare <- expand.grid(
  x1 = seq(0, 70, by = 0.25),
  x2 = seq(0, 50, by = 0.25)
)

grid_compare$Original <- mapply(is_feasible, grid_compare$x1, grid_compare$x2)
grid_compare$Modified <- mapply(is_feasible_modified, grid_compare$x1, grid_compare$x2)

grid_original <- subset(grid_compare, Original)
grid_modified <- subset(grid_compare, Modified)

x_mod <- seq(0, 70, length.out = 400)
line1_old <- 100 - 2 * x_mod
line1_new <- 120 - 2 * x_mod
line2_mod <- (80 - x_mod) / 2

p6 <- ggplot() +
  geom_point(
    data = grid_modified,
    aes(x = x1, y = x2),
    alpha = 0.08,
    size = 0.5
  ) +
  geom_point(
    data = grid_original,
    aes(x = x1, y = x2),
    alpha = 0.15,
    size = 0.5
  ) +
  geom_line(
    data = data.frame(x = x_mod, y = line1_old),
    aes(x = x, y = y),
    linewidth = 1,
    linetype = "dashed"
  ) +
  geom_line(
    data = data.frame(x = x_mod, y = line1_new),
    aes(x = x, y = y),
    linewidth = 1
  ) +
  geom_line(
    data = data.frame(x = x_mod, y = line2_mod),
    aes(x = x, y = y),
    linewidth = 1
  ) +
  labs(
    title = "Sensitivity: Original vs Relaxed Constraint",
    subtitle = "Dashed line = original constraint, solid line = relaxed constraint",
    x = "x1 (Product A)",
    y = "x2 (Product B)"
  ) +
  coord_cartesian(xlim = c(0, 70), ylim = c(0, 50)) +
  theme_minimal(base_size = 13)

print(p6)

############################
# 5. SECOND EXAMPLE WITH NON-ZERO SLACK
############################

cat("\n==============================\n")
cat("SECOND EXAMPLE: NON-ZERO SLACK\n")
cat("==============================\n")

# This second example is included because in the first example
# both constraints are binding, so slack is zero.
#
# Maximise:
# z = 8x1 + 11x2
#
# Constraints:
# 3x1 + 5x2 <= 150   (budget)
# 2x1 + x2 <= 100    (team hours)
# x1, x2 >= 0

objective2 <- c(8, 11)

constraint_matrix2 <- matrix(c(
  3, 5,
  2, 1
), nrow = 2, byrow = TRUE)

constraint_rhs2 <- c(150, 100)

lp_result2 <- lpSolve::lp(
  direction = "max",
  objective.in = objective2,
  const.mat = constraint_matrix2,
  const.dir = c("<=", "<="),
  const.rhs = constraint_rhs2,
  compute.sens = TRUE
)

cat("\nOptimal objective value (example 2):\n")
print(lp_result2$objval)

cat("\nOptimal decision variables (example 2):\n")
print(lp_result2$solution)

used_resources2 <- as.numeric(constraint_matrix2 %*% lp_result2$solution)
slack2 <- constraint_rhs2 - used_resources2
dual_prices2 <- lp_result2$duals[1:2]

results_table2 <- data.frame(
  Constraint = c("Budget", "Team hours"),
  Capacity = constraint_rhs2,
  Used = used_resources2,
  Slack = slack2,
  ShadowPrice = dual_prices2
)

cat("\nConstraint usage, slack, and shadow prices (example 2):\n")
print(results_table2)

# Plot 7: Slack in the second example
p7 <- ggplot(results_table2, aes(x = Constraint, y = Slack)) +
  geom_col() +
  geom_text(aes(label = round(Slack, 2)), vjust = -0.4, size = 4) +
  labs(
    title = "Slack in the Second Example",
    subtitle = "This example is useful because one constraint can have visible slack",
    x = "Constraint",
    y = "Slack"
  ) +
  theme_minimal(base_size = 13)

print(p7)

# Plot 8: Resource usage in the second example
resource_plot_data2 <- rbind(
  data.frame(Constraint = results_table2$Constraint, Type = "Capacity", Value = results_table2$Capacity),
  data.frame(Constraint = results_table2$Constraint, Type = "Used", Value = results_table2$Used)
)

p8 <- ggplot(resource_plot_data2, aes(x = Constraint, y = Value, fill = Type)) +
  geom_col(position = "dodge") +
  geom_text(
    aes(label = round(Value, 1)),
    position = position_dodge(width = 0.9),
    vjust = -0.4,
    size = 4
  ) +
  labs(
    title = "Resource Usage vs Capacity (Second Example)",
    x = "Constraint",
    y = "Units"
  ) +
  theme_minimal(base_size = 13)

print(p8)

############################
# 6. SUMMARY
############################

cat("\n==============================\n")
cat("SUMMARY\n")
cat("==============================\n")

cat("
1. A decision problem becomes a mathematical model.
2. The feasible region contains all valid decisions.
3. In a two-variable LP, the optimum occurs at a corner point.
4. Binding constraints are fully used.
5. Slack measures unused resources.
6. Shadow prices measure the value of an extra unit of a resource.
7. Plots help us interpret the model, not only solve it.
\n")