# ============================================================
# Programmer: Daniel Hriscu
# Student ID: 2127929
# Title: Alien Colony Resource Management Simulation
# Module: ACCA7022 - Advanced Statistical Programming
# ============================================================

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

# Set seed for reproducibility so results can be independently verified
set.seed(2127929)

# Define constants to avoid "magic numbers" in the code.
# Centralising these parameters makes the simulation easy to
# reconfigure and reduces the risk of inconsistent hard-coded values.
base_output  <- 50
num_sols     <- 30
num_stations <- 6
num_teams    <- 4

# INPUT VALIDATION
# Verify simulation parameters are within expected bounds before
# proceeding, to prevent silent logical errors at runtime.
if (base_output <= 0) stop("base_output must be positive")
if (num_sols < 1 || num_teams < 1 || num_stations < 1) {
  stop("Simulation dimensions must be positive integers")
}

# Calculate total expected rows (30 sols * 6 stations * 4 teams = 720)
total_rows <- num_sols * num_stations * num_teams

# Pre-allocate the dataframe with typed vectors.
# This avoids the O(n^2) penalty of growing a dataframe with rbind(),
# which forces R to copy the entire object on every append [1].
production_data <- data.frame(
  Sol     = integer(total_rows),
  Station = integer(total_rows),
  Team    = integer(total_rows),
  Output  = numeric(total_rows)
)

# Populate data using nested loops with a row_index counter.
# The loop order (sol -> station -> team) creates a deterministic
# row layout that later tasks exploit for direct index arithmetic.
row_index <- 1
for (sol in 1:num_sols) {
  for (station in 1:num_stations) {
    for (team in 1:num_teams) {

      # Generate random noise: discrete uniform integer in [-5, +5]
      random_noise <- sample(-5:5, 1)

      # Calculate daily production output using the specified formula
      output_val <- base_output + (station * 3) + (team * 2) + random_noise

      # Assign values directly to the pre-allocated dataframe via row index
      production_data$Sol[row_index]     <- sol
      production_data$Station[row_index] <- station
      production_data$Team[row_index]    <- team
      production_data$Output[row_index]  <- output_val

      # Increment the index for the next iteration
      row_index <- row_index + 1
    }
  }
}

# Display the first 10 rows of the dataframe
cat("\n--- TASK 1A: First 10 Rows of Production Data ---\n\n")
print(head(production_data, 10), row.names = FALSE)

# ============================================================
# TESTING TASK 1A
# ============================================================
cat("\n--- TASK 1A TESTING ---\n")

# TEST 1: Point Verification
cat("Row 1 Output:", production_data$Output[1], "\n")
cat("Expected base (no noise): 55 | Valid range: [50, 60]\n")
cat("Within bounds?",
    production_data$Output[1] >= 50 & production_data$Output[1] <= 60, "\n\n")

# TEST 2: Structural Validation (loop-based unique counting)
unique_sols <- 0
seen_sols <- integer(0)
unique_stations <- 0
seen_stations <- integer(0)
unique_teams <- 0
seen_teams <- integer(0)

for (i in 1:nrow(production_data)) {
  is_new_sol <- TRUE
  if (unique_sols > 0) {
    for (j in 1:unique_sols) {
      if (seen_sols[j] == production_data$Sol[i]) is_new_sol <- FALSE
    }
  }
  if (is_new_sol) {
    unique_sols <- unique_sols + 1
    seen_sols[unique_sols] <- production_data$Sol[i]
  }

  is_new_station <- TRUE
  if (unique_stations > 0) {
    for (j in 1:unique_stations) {
      if (seen_stations[j] == production_data$Station[i]) is_new_station <- FALSE
    }
  }
  if (is_new_station) {
    unique_stations <- unique_stations + 1
    seen_stations[unique_stations] <- production_data$Station[i]
  }

  is_new_team <- TRUE
  if (unique_teams > 0) {
    for (j in 1:unique_teams) {
      if (seen_teams[j] == production_data$Team[i]) is_new_team <- FALSE
    }
  }
  if (is_new_team) {
    unique_teams <- unique_teams + 1
    seen_teams[unique_teams] <- production_data$Team[i]
  }
}

cat("Total Rows:", nrow(production_data), "(Expected: 720)\n")
cat("Unique Sols:", unique_sols, "(Expected: 30)\n")
cat("Unique Stations:", unique_stations, "(Expected: 6)\n")
cat("Unique Teams:", unique_teams, "(Expected: 4)\n")
cat("Rows per Sol-Station combo:", total_rows / (num_sols * num_stations),
    "(Expected: 4)\n\n")

# TEST 3: Range Validation (loop-based min/max)
theoretical_min <- base_output + (1 * 3) + (1 * 2) + (-5)
theoretical_max <- base_output + (num_stations * 3) + (num_teams * 2) + 5

observed_min <- production_data$Output[1]
observed_max <- production_data$Output[1]
for (i in 2:nrow(production_data)) {
  if (production_data$Output[i] < observed_min) observed_min <- production_data$Output[i]
  if (production_data$Output[i] > observed_max) observed_max <- production_data$Output[i]
}

cat("Min Output:", observed_min, "(Theoretical min:", theoretical_min, ")\n")
cat("Max Output:", observed_max, "(Theoretical max:", theoretical_max, ")\n")

all_in_range <- TRUE
for (i in 1:nrow(production_data)) {
  if (production_data$Output[i] < theoretical_min ||
      production_data$Output[i] > theoretical_max) {
    all_in_range <- FALSE
  }
}
cat("All outputs within theoretical bounds?", all_in_range, "\n")


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

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

# SINGLE-PASS ACCUMULATION: O(n) time complexity.
# Rather than a nested approach (outer loop over stations, inner loop
# scanning all rows with an if-condition, giving O(n*k) comparisons),
# we iterate exactly once through 720 rows. Each row's Station ID is
# used as a direct index into the accumulator vector, eliminating all
# conditional branching [2].
for (i in 1:nrow(production_data)) {
  current_station <- production_data$Station[i]
  current_output  <- production_data$Output[i]

  station_totals$TotalOutput[current_station] <-
    station_totals$TotalOutput[current_station] + current_output
}

cat("\n--- TASK 1B: Total Production per Station ---\n\n")
print(station_totals, row.names = FALSE)

# ============================================================
# TESTING TASK 1B
# ============================================================
cat("\n--- TASK 1B TESTING ---\n")

# TEST 1: Boundary Verification (Station 1 and Station 6)
manual_sum_station1 <- 0
manual_sum_station6 <- 0
for (i in 1:nrow(production_data)) {
  if (production_data$Station[i] == 1) {
    manual_sum_station1 <- manual_sum_station1 + production_data$Output[i]
  }
  if (production_data$Station[i] == 6) {
    manual_sum_station6 <- manual_sum_station6 + production_data$Output[i]
  }
}

cat("Station 1 Expected (manual loop):", manual_sum_station1, "\n")
cat("Station 1 Actual   (accumulator):", station_totals$TotalOutput[1], "\n")
cat("Station 1 TEST PASSED?", manual_sum_station1 == station_totals$TotalOutput[1], "\n\n")

cat("Station 6 Expected (manual loop):", manual_sum_station6, "\n")
cat("Station 6 Actual   (accumulator):", station_totals$TotalOutput[6], "\n")
cat("Station 6 TEST PASSED?", manual_sum_station6 == station_totals$TotalOutput[6], "\n\n")

# TEST 2: Cross-Validation (Grand Total)
grand_total_from_stations <- 0
for (s in 1:num_stations) {
  grand_total_from_stations <- grand_total_from_stations + station_totals$TotalOutput[s]
}
grand_total_from_data <- 0
for (i in 1:nrow(production_data)) {
  grand_total_from_data <- grand_total_from_data + production_data$Output[i]
}
cat("Grand Total (station_totals):", grand_total_from_stations, "\n")
cat("Grand Total (production_data):", grand_total_from_data, "\n")
cat("Grand Total TEST PASSED?", grand_total_from_stations == grand_total_from_data, "\n")


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

anomalies <- production_data
anomalies$Anomaly <- FALSE

# CHUNKED TWO-PASS APPROACH: O(n) total complexity.
# Because the data is deterministically ordered (4 teams per station
# per sol), we process sequential chunks of 4 rows using a pointer.
# Pass 1: aggregate the sum to compute the daily station average.
# Pass 2: evaluate each team against the 10% deviation threshold.
# The mean cannot be known until all 4 values are read, so the
# two-pass design is a sequential dependency, not inefficiency [3].
current_idx <- 1

for (sol in 1:num_sols) {
  for (station in 1:num_stations) {

    # PASS 1: Sum the 4 team outputs to compute the daily average
    daily_sum <- 0
    for (team_offset in 0:(num_teams - 1)) {
      daily_sum <- daily_sum + anomalies$Output[current_idx + team_offset]
    }
    daily_avg <- daily_sum / num_teams

    # PASS 2: Flag teams deviating by more than 10% from the daily average
    for (team_offset in 0:(num_teams - 1)) {
      row_to_check <- current_idx + team_offset
      team_output  <- anomalies$Output[row_to_check]

      if (abs(team_output - daily_avg) > (0.10 * daily_avg)) {
        anomalies$Anomaly[row_to_check] <- TRUE
      }
    }

    current_idx <- current_idx + num_teams
  }
}

cat("\n--- TASK 1C: First 10 Anomalies Detected ---\n\n")

first_10_anomalies <- data.frame(
  Sol = integer(10), Station = integer(10), Team = integer(10),
  Output = numeric(10), Anomaly = logical(10)
)

found_count <- 0
for (i in 1:nrow(anomalies)) {
  if (anomalies$Anomaly[i] == TRUE) {
    found_count <- found_count + 1
    first_10_anomalies$Sol[found_count]     <- anomalies$Sol[i]
    first_10_anomalies$Station[found_count] <- anomalies$Station[i]
    first_10_anomalies$Team[found_count]    <- anomalies$Team[i]
    first_10_anomalies$Output[found_count]  <- anomalies$Output[i]
    first_10_anomalies$Anomaly[found_count] <- anomalies$Anomaly[i]
  }
  if (found_count == 10) break
}

if (found_count < 10) {
  print(first_10_anomalies[1:found_count, ], row.names = FALSE)
} else {
  print(first_10_anomalies, row.names = FALSE)
}

# ============================================================
# TESTING TASK 1C
# ============================================================
cat("\n--- TASK 1C TESTING ---\n")

# TEST 1: Point Verification (Sol 1, Station 1)
test_sum <- 0
for (t in 1:4) test_sum <- test_sum + production_data$Output[t]
manual_avg <- test_sum / 4
threshold  <- manual_avg * 0.10

cat("Sol 1, Station 1 Outputs:",
    production_data$Output[1], production_data$Output[2],
    production_data$Output[3], production_data$Output[4], "\n")
cat("Manual Daily Average:", manual_avg, "\n")
cat("10% Threshold (Allowed Deviation):", threshold, "\n")

for (t in 1:4) {
  dev  <- abs(production_data$Output[t] - manual_avg)
  flag <- dev > threshold
  cat("  Team", t, "| Output:", production_data$Output[t],
      "| Deviation:", dev, "| Anomaly:", flag, "\n")
}

# TEST 2: Negative Case Verification
cat("\nRow 2 Anomaly flag:", anomalies$Anomaly[2],
    "(Expected: FALSE based on manual check above)\n")

# TEST 3: Statistical Reasonableness
total_anomalies <- 0
for (i in 1:nrow(anomalies)) {
  if (anomalies$Anomaly[i]) total_anomalies <- total_anomalies + 1
}
anomaly_rate <- round(total_anomalies / nrow(anomalies) * 100, 1)
cat("Total anomalies flagged:", total_anomalies, "out of", nrow(anomalies), "\n")
cat("Anomaly rate:", anomaly_rate, "%\n")


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

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

# DIRECT INDEX ARITHMETIC: For each station, iterate through all 30 sols.
# Instead of scanning the full dataframe, the exact starting row for any
# Sol/Station combination is computed using the deterministic data layout:
#   Index = ((sol-1) * num_stations * num_teams) + ((station-1) * num_teams) + 1
# This allows O(1) random access per chunk, giving O(s*d) total [3].
for (target_station in 1:num_stations) {

  max_output <- -Inf
  best_sol   <- 0

  for (target_sol in 1:num_sols) {
    current_sol_total <- 0

    start_idx <- ((target_sol - 1) * num_stations * num_teams) +
                 ((target_station - 1) * num_teams) + 1

    for (team_offset in 0:(num_teams - 1)) {
      current_sol_total <- current_sol_total +
        production_data$Output[start_idx + team_offset]
    }

    if (current_sol_total > max_output) {
      max_output <- current_sol_total
      best_sol   <- target_sol
    }
  }

  peak_days$Sol[target_station]         <- best_sol
  peak_days$TotalOutput[target_station] <- max_output
}

cat("\n--- TASK 1D: Peak Productivity Days ---\n\n")
print(peak_days, row.names = FALSE)

# ============================================================
# TESTING TASK 1D
# ============================================================
cat("\n--- TASK 1D TESTING ---\n")

# TEST 1: Independent Verification for Station 1 (loop-based)
best_sol_test1   <- 0
best_total_test1 <- -Inf

for (s in 1:num_sols) {
  sol_total <- 0
  idx <- ((s - 1) * num_stations * num_teams) + ((1 - 1) * num_teams) + 1
  for (t in 0:(num_teams - 1)) {
    sol_total <- sol_total + production_data$Output[idx + t]
  }
  if (sol_total > best_total_test1) {
    best_total_test1 <- sol_total
    best_sol_test1   <- s
  }
}

cat("Station 1 Expected Peak Sol:", best_sol_test1, "with Output:", best_total_test1, "\n")
cat("Station 1 Loop Found Sol:   ", peak_days$Sol[1], "with Output:", peak_days$TotalOutput[1], "\n")
cat("Station 1 TEST PASSED?",
    best_sol_test1 == peak_days$Sol[1] && best_total_test1 == peak_days$TotalOutput[1], "\n\n")

# TEST 2: Cross-Validate Station 3
best_sol_test3   <- 0
best_total_test3 <- -Inf

for (s in 1:num_sols) {
  sol_total <- 0
  idx <- ((s - 1) * num_stations * num_teams) + ((3 - 1) * num_teams) + 1
  for (t in 0:(num_teams - 1)) {
    sol_total <- sol_total + production_data$Output[idx + t]
  }
  if (sol_total > best_total_test3) {
    best_total_test3 <- sol_total
    best_sol_test3   <- s
  }
}

cat("Station 3 Expected Peak Sol:", best_sol_test3, "with Output:", best_total_test3, "\n")
cat("Station 3 Loop Found Sol:   ", peak_days$Sol[3], "with Output:", peak_days$TotalOutput[3], "\n")
cat("Station 3 TEST PASSED?",
    best_sol_test3 == peak_days$Sol[3] && best_total_test3 == peak_days$TotalOutput[3], "\n\n")

# TEST 3: Sanity Check
all_positive <- TRUE
for (s in 1:num_stations) {
  if (peak_days$TotalOutput[s] <= 0) all_positive <- FALSE
}
cat("All peak totals > 0?", all_positive, "\n")


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

# Pre-allocate the new column as a numeric vector to avoid
# memory fragmentation from incremental column growth [1].
production_data$EfficiencyScore <- numeric(nrow(production_data))

# EDGE CASE ANALYSIS: The denominator (Station + Team) has a theoretical
# range of [2, 10] (min: Station=1+Team=1, max: Station=6+Team=4).
# Division by zero is therefore mathematically impossible within this
# simulation's parameter space. EfficiencyScore is bounded by [5.0, 40.5].
for (i in 1:nrow(production_data)) {
  current_output  <- production_data$Output[i]
  current_station <- production_data$Station[i]
  current_team    <- production_data$Team[i]

  production_data$EfficiencyScore[i] <- current_output / (current_station + current_team)
}

cat("\n--- TASK 1E: First 10 Rows with Efficiency Score ---\n\n")
print(head(production_data, 10), row.names = FALSE)

# ============================================================
# TESTING TASK 1E
# ============================================================
cat("\n--- TASK 1E TESTING ---\n")

# TEST 1: Point Verification (Row 1)
row1_output  <- production_data$Output[1]
row1_station <- production_data$Station[1]
row1_team    <- production_data$Team[1]
manual_efficiency_row1 <- row1_output / (row1_station + row1_team)
loop_efficiency_row1   <- production_data$EfficiencyScore[1]
cat("Row 1 Data: Output =", row1_output,
    ", Station =", row1_station, ", Team =", row1_team, "\n")
cat("Expected Efficiency:", manual_efficiency_row1, "\n")
cat("Loop Calculated:    ", loop_efficiency_row1, "\n")
cat("Row 1 TEST PASSED?", manual_efficiency_row1 == loop_efficiency_row1, "\n\n")

# TEST 2: Range Validation
eff_min <- production_data$EfficiencyScore[1]
eff_max <- production_data$EfficiencyScore[1]
for (i in 2:nrow(production_data)) {
  if (production_data$EfficiencyScore[i] < eff_min) eff_min <- production_data$EfficiencyScore[i]
  if (production_data$EfficiencyScore[i] > eff_max) eff_max <- production_data$EfficiencyScore[i]
}
cat("Min EfficiencyScore:", eff_min, "(Theoretical min: 5.0)\n")
cat("Max EfficiencyScore:", eff_max, "(Theoretical max: 40.5)\n")
cat("All within bounds?", eff_min >= 5.0 & eff_max <= 40.5, "\n")