# ==============================================================================
# ACCA7022: Alien Colony Resource Management 
# ==============================================================================

# --- Setup Variables ---
total_sols <- 30
total_stations <- 6
total_teams <- 4
base_output <- 50

total_records <- total_sols * total_stations * total_teams
sol_vec <- numeric(total_records)
station_vec <- numeric(total_records)
team_vec <- numeric(total_records)
output_vec <- numeric(total_records)

row_index <- 1

# ==============================================================================
# TASK 1A: Generate Production Data
# ==============================================================================

for (s in 1:total_sols) {
  for (st in 1:total_stations) {
    for (t in 1:total_teams) {
      
      sol_vec[row_index] <- s
      station_vec[row_index] <- st
      team_vec[row_index] <- t
      
      # Output = BaseOutput + (Station * 3) + (Team * 2) + RandomNoise
      random_noise <- sample(-5:5, 1)
      daily_output <- base_output + (st * 3) + (t * 2) + random_noise
      output_vec[row_index] <- daily_output
      
      row_index <- row_index + 1
    }
  }
}

production_data <- data.frame(
  Sol = sol_vec, 
  Station = station_vec, 
  Team = team_vec, 
  Output = output_vec
)

# --- Console Output 1A ---
cat("\n======================================================\n")
cat(" TASK 1A: GENERATED PRODUCTION DATA (FIRST 10 ROWS)\n")
cat("======================================================\n")
print(head(production_data, 10))
cat("------------------------------------------------------\n")

# --- Graphs for 1A ---
# Graph 1: Distribution of all production outputs
hist(production_data$Output, breaks = 15, col = "skyblue", border = "white",
     main = "Task 1A: Distribution of Daily Outputs",
     xlab = "Output Value", ylab = "Frequency")

# Graph 2: Output spread broken down by Worker Teams
boxplot(Output ~ Team, data = production_data, col = c("gold", "tomato", "lightgreen", "plum"),
        main = "Task 1A: Output Spread by Worker Team",
        xlab = "Team ID", ylab = "Daily Output")


# ==============================================================================
# TASK 1B: Total Production per Station
# ==============================================================================

station_totals <- data.frame(
  Station = 1:total_stations, 
  TotalOutput = numeric(total_stations)
)

for (i in 1:nrow(production_data)) {
  current_station <- production_data$Station[i]
  current_output <- production_data$Output[i]
  
  # Accumulate totals for the matching station
  station_totals$TotalOutput[station_totals$Station == current_station] <- 
    station_totals$TotalOutput[station_totals$Station == current_station] + current_output
}

# --- Console Output 1B ---
cat("\n======================================================\n")
cat(" TASK 1B: TOTAL PRODUCTION PER STATION\n")
cat("======================================================\n")
print(station_totals, row.names = FALSE)
cat("------------------------------------------------------\n")

# --- Graphs for 1B ---
# Graph 3: Bar plot of total outputs per station
barplot(station_totals$TotalOutput, names.arg = paste("Station", station_totals$Station),
        col = "steelblue", main = "Task 1B: Total Output by Station",
        xlab = "Station ID", ylab = "Total Output (Units)")

# Graph 4: Pie chart showing the proportion of total output
pie(station_totals$TotalOutput, labels = paste("St", station_totals$Station), 
    main = "Task 1B: Share of Total Output by Station", col = rainbow(6))


# ==============================================================================
# TASK 1C: Detect Production Anomalies
# ==============================================================================

daily_sums <- matrix(0, nrow = total_sols, ncol = total_stations)

for (i in 1:nrow(production_data)) {
  s <- production_data$Sol[i]
  st <- production_data$Station[i]
  daily_sums[s, st] <- daily_sums[s, st] + production_data$Output[i]
}

anomaly_vec <- logical(nrow(production_data))
color_vec <- character(nrow(production_data)) 

for (i in 1:nrow(production_data)) {
  s <- production_data$Sol[i]
  st <- production_data$Station[i]
  out <- production_data$Output[i]
  
  daily_avg <- daily_sums[s, st] / 4
  
  # 10% bounds
  lower_bound <- daily_avg * 0.90
  upper_bound <- daily_avg * 1.10
  
  if (out < lower_bound || out > upper_bound) {
    anomaly_vec[i] <- TRUE
    color_vec[i] <- "red" # Highlight anomalies in red
  } else {
    anomaly_vec[i] <- FALSE
    color_vec[i] <- "darkgray" # Normal data in gray
  }
}

production_data$Anomaly <- anomaly_vec
anomalies <- production_data[production_data$Anomaly == TRUE, c("Sol", "Station", "Team", "Output", "Anomaly")]

# --- Console Output 1C ---
cat("\n======================================================\n")
cat(" TASK 1C: DETECTED PRODUCTION ANOMALIES (FIRST 10)\n")
cat("======================================================\n")
print(head(anomalies, 10))
cat("------------------------------------------------------\n")

# --- Graphs for 1C ---
# Graph 5: Scatter plot highlighting anomalous outputs in red
plot(production_data$Output, col = color_vec, pch = 16,
     main = "Task 1C: Production Outputs (Anomalies in Red)",
     xlab = "Observation Index", ylab = "Daily Output")


# ==============================================================================
# TASK 1D: Identify Peak Productivity Days
# ==============================================================================

peak_days <- data.frame(Station = integer(), Sol = integer(), TotalOutput = numeric())

for (st in 1:total_stations) {
  max_output <- -1
  peak_sol <- -1
  
  for (s in 1:total_sols) {
    current_day_total <- daily_sums[s, st]
    
    if (current_day_total > max_output) {
      max_output <- current_day_total
      peak_sol <- s
    }
  }
  
  peak_days[nrow(peak_days) + 1, ] <- c(st, peak_sol, max_output)
}

# --- Console Output 1D ---
cat("\n======================================================\n")
cat(" TASK 1D: PEAK PRODUCTIVITY DAYS PER STATION\n")
cat("======================================================\n")
print(peak_days, row.names = FALSE)
cat("------------------------------------------------------\n")

# --- Graphs for 1D ---
# Graph 6: Bar plot comparing the peak day outputs across stations
barplot(peak_days$TotalOutput, names.arg = paste("Station", peak_days$Station),
        col = "darkorange", main = "Task 1D: Peak Daily Output per Station",
        xlab = "Station ID", ylab = "Peak Day Total Output")


# ==============================================================================
# TASK 1E: Calculate Efficiency Score
# ==============================================================================

efficiency_vec <- numeric(nrow(production_data))

for (i in 1:nrow(production_data)) {
  out <- production_data$Output[i]
  st <- production_data$Station[i]
  t <- production_data$Team[i]
  
  efficiency_vec[i] <- out / (st + t)
}

production_data$EfficiencyScore <- efficiency_vec

# --- Console Output 1E ---
cat("\n======================================================\n")
cat(" TASK 1E: EFFICIENCY SCORES (FIRST 10 ROWS)\n")
cat("======================================================\n")
print(head(production_data, 10))
cat("======================================================\n")

# --- Graphs for 1E ---
# Graph 7: Histogram of efficiency scores
hist(production_data$EfficiencyScore, breaks = 20, col = "mediumpurple", border = "white",
     main = "Task 1E: Distribution of Efficiency Scores",
     xlab = "Efficiency Score", ylab = "Frequency")

# Graph 8: Scatter plot showing Output vs Efficiency Score
plot(production_data$Output, production_data$EfficiencyScore,
     main = "Task 1E: Output vs Efficiency Score",
     xlab = "Daily Output", ylab = "Efficiency Score",
     pch = 19, col = rgb(0.2, 0.6, 0.4, 0.5))