
# Smart City Utilities Operations Simulation


set.seed(42)  # Set seed for reproducibility


# TASK A: Generate Production Data

# Define simulation parameters
n_days      <- 30
n_utilities <- 6
n_teams     <- 5
base_output <- 40

# Pre-calculate total number of rows: 30 days x 6 utilities x 5 teams
total_rows <- n_days * n_utilities * n_teams

# Initialise empty vectors to hold each column's data
# Pre-allocating is more efficient than growing a data frame row by row
day_col     <- integer(total_rows)
utility_col <- integer(total_rows)
team_col    <- integer(total_rows)
output_col  <- numeric(total_rows)

# Row index counter — incremented inside the innermost loop
row_index <- 1

# Nested loops
for (day in 1:n_days) {
  for (utility in 1:n_utilities) {
    for (team in 1:n_teams) {
      
      # Generate random noise: one integer uniformly drawn from -5 to +5
      random_noise <- sample(-5:5, size = 1)
      
      # Calculate daily output using the specified formula:
      # Output = BaseOutput + (Utility x 3) + (Team x 2) + RandomNoise
      output_value <- base_output + (utility * 3) + (team * 2) + random_noise
      
      # Store computed values in their respective column vectors
      day_col[row_index]     <- day
      utility_col[row_index] <- utility
      team_col[row_index]    <- team
      output_col[row_index]  <- output_value
      
      # Advance the row index before the next iteration
      row_index <- row_index + 1
    }
  }
}

# Assemble the four vectors into a single data frame
production_data <- data.frame(
  Day     = day_col,
  Utility = utility_col,
  Team    = team_col,
  Output  = output_col
)

# Display the first 10 rows of the production data
cat("=== Task A: First 10 Rows of production_data ===\n")
print(head(production_data, 10))
cat("\n")



# TASK B: Total Production per Utility


# Initialise a vector to accumulate the total output for each utility
utility_total_output <- numeric(n_utilities)

# Loop over every row of production_data to sum output by utility
for (i in 1:nrow(production_data)) {
  
  # Read the utility ID for this row
  u <- production_data$Utility[i]
  
  # Add this row's output to the running total for utility u
  utility_total_output[u] <- utility_total_output[u] + production_data$Output[i]
}

# Build the Utility_totals data frame from the accumulated sums
Utility_totals <- data.frame(
  Utility     = 1:n_utilities,
  TotalOutput = utility_total_output
)

# Display total output per utility
cat("=== Task B: Total Production per Utility ===\n")
print(Utility_totals)
cat("\n")



# TASK C: Detect Production Anomalies

# Build a [day x utility] matrix of daily total outputs
# Using the matrix avoids forbidden aggregation functions
daily_utility_sum <- matrix(0, nrow = n_days, ncol = n_utilities)

for (i in 1:nrow(production_data)) {
  d <- production_data$Day[i]
  u <- production_data$Utility[i]
  daily_utility_sum[d, u] <- daily_utility_sum[d, u] + production_data$Output[i]
}

# Compute daily average per utility (sum divided by 4, as specified)
daily_utility_avg <- matrix(0, nrow = n_days, ncol = n_utilities)

for (d in 1:n_days) {
  for (u in 1:n_utilities) {
    daily_utility_avg[d, u] <- daily_utility_sum[d, u] / 4
  }
}

# Populate the anomalies data frame row by row
anom_day     <- integer(total_rows)
anom_utility <- integer(total_rows)
anom_team    <- integer(total_rows)
anom_output  <- numeric(total_rows)
anom_flag    <- logical(total_rows)

for (i in 1:nrow(production_data)) {
  d      <- production_data$Day[i]
  u      <- production_data$Utility[i]
  t      <- production_data$Team[i]
  output <- production_data$Output[i]
  avg    <- daily_utility_avg[d, u]
  
  # Percentage deviation of this team's output from the daily utility average
  pct_deviation <- abs(output - avg) / avg
  
  # Anomaly is TRUE when deviation exceeds 10%
  is_anomaly <- pct_deviation > 0.10
  
  # Store all values for this row
  anom_day[i]     <- d
  anom_utility[i] <- u
  anom_team[i]    <- t
  anom_output[i]  <- output
  anom_flag[i]    <- is_anomaly
}

# Assemble the anomalies data frame
anomalies <- data.frame(
  Day     = anom_day,
  Utility = anom_utility,
  Team    = anom_team,
  Output  = anom_output,
  Anomaly = anom_flag
)

# Filter to rows where Anomaly == TRUE and display the first 10
anomaly_rows <- anomalies[anomalies$Anomaly == TRUE, ]

cat("=== Task C: First 10 Detected Anomalies ===\n")
print(head(anomaly_rows, 10))


# TASK D: Identify Peak Productivity Days

# Initialise vectors to store the peak day and peak total output per utility
peak_day_vec    <- integer(n_utilities)
peak_output_vec <- numeric(n_utilities)

for (u in 1:n_utilities) {
  
  # Start with day 1 as the provisional best for this utility
  best_day    <- 1
  best_output <- daily_utility_sum[1, u]
  
  # Scan remaining days, updating the best whenever a higher total is found
  for (d in 2:n_days) {
    if (daily_utility_sum[d, u] > best_output) {
      best_output <- daily_utility_sum[d, u]
      best_day    <- d
    }
  }
  
  peak_day_vec[u]    <- best_day
  peak_output_vec[u] <- best_output
}

# Build the peak_days data frame
peak_days <- data.frame(
  Utility     = 1:n_utilities,
  Day         = peak_day_vec,
  TotalOutput = peak_output_vec
)

cat("=== Task D: Peak Productivity Day per Utility ===\n")
print(peak_days)
cat("\n")



# TASK E: Calculate Efficiency Score

# Pre-allocate the efficiency score vector
efficiency_scores <- numeric(total_rows)

for (i in 1:nrow(production_data)) {
  u <- production_data$Utility[i]
  t <- production_data$Team[i]
  
  # Efficiency score: actual output divided by (Utility + Team)
  efficiency_scores[i] <- production_data$Output[i] / (u + t)
}

# Append the new column to production_data
production_data$EfficiencyScore <- efficiency_scores

# Display the first 10 rows of the updated data frame
cat("=== Task E: First 10 Rows with EfficiencyScore ===\n")
print(head(production_data, 10))
cat("\n")


