### Demographic Model ###



# -------------------------------
# Load required libraries
# -------------------------------
library(ggplot2)   # Plotting
library(writexl)   # Exporting Excel files
library(dplyr)     # Data manipulation
library(tidyr)     # Data reshaping


# -------------------------------
# Simulation Parameters
# -------------------------------
set.seed(1234)  # Ensures reproducibility for simulations

n_each_model <- 1000000  # Number of simulation runs per model type

# Time range of possible settlement duration (in years)
duration_range <- c(123, 123)

# Observed archaeological constraints from excavation evidence
observed_min <- 123    # Minimum number of estimated households
observed_max <- 123    # Maximum number of estimated households
duration_max <- 123    # Maximum settlement duration allowed

# Ranges for simulation input parameters
initial_pop_range <- c(123, 123)              # Starting population size
carrying_capacity_range <- c(123, 123)      # Maximum sustainable population
growth_rate_range_pos <- c(123, 123)           # Positive growth rate range (e.g., expansion phase)
growth_rate_range_neg <- c(123, 123)           # Negative growth rate range (e.g., abandonment)
household_size_range <- c(123, 123)               # Number of people per household
house_lifetime_range <- c(123, 123)             # Average lifespan of a house in years


# -------------------------------
# Functions
# -------------------------------

# Simulates logistic growth (used for both models)
logistic_growth <- function(P0, K, r, T) {
  t <- 0:(T - 1)
  K / (1 + ((K - P0) / P0) * exp(-r * t))
}

# Estimates the effective duration (i.e., time when population is at least 1% of its peak)
estimate_duration <- function(P) {
  max_P <- max(P)
  for (i in seq_along(P)) {
    if (P[i] >= 0.99 * max_P || P[i] <= 0.01 * max_P) return(i)
  }
  return(length(P))
}


# -------------------------------
# Initialize Results Storage
# -------------------------------
results_ram <- data.frame()  # Stores valid simulations from Rapid Abandonment Model
results_gam <- data.frame()  # Stores valid simulations from Gradual Abandonment Model
curves_ram <- list()         # Stores population trajectories for RAM
curves_gam <- list()         # Stores population trajectories for GAM


# -------------------------------
# Rapid Abandonment Model 
# -------------------------------

for (i in 1:n_each_model) {
  T <- sample(duration_range[1]:duration_range[2], 1)
  t <- 0:(T - 1)
  
  # Randomly sample parameters
  P0 <- runif(1, initial_pop_range[1], initial_pop_range[2])
  K <- runif(1, max(P0 + 1, carrying_capacity_range[1]), carrying_capacity_range[2])
  r <- runif(1, growth_rate_range_pos[1], growth_rate_range_pos[2])
  S <- runif(1, household_size_range[1], household_size_range[2])
  L <- rnorm(1, mean = mean(house_lifetime_range), sd = 5)
  
  # Skip invalid inputs
  if (L < house_lifetime_range[1] || L > house_lifetime_range[2] || S <= 0 || L <= 0) next

  # Simulate population curve
  P <- logistic_growth(P0, K, r, T)
  max_pop <- max(P)
  if (abs(K - max_pop) >= 10) next
  
  # Calculate number of households over the site lifespan
  steps <- floor(T / L)
  house_counts <- sapply(0:(steps - 1), function(i) {
    index <- i * floor(L) + 1
    if (index <= length(P)) return(P[index] / S) else return(0)
  })
  H <- sum(unlist(house_counts))
  duration <- estimate_duration(P)
  
  # Enforce constraints: initial houses should not exceed 50
  if (P[1] / S > 50) next

  # Final constraint: match observed archaeological household range
  if (H >= observed_min && H <= observed_max && duration <= duration_max) {
   # Write results
   results_ram <- rbind(results_ram, data.frame(
      Initial_Population = P[1],
      Final_Population = P[length(P)],
      Max_Population = max_pop,
      Growth_Rate = r,
      Duration = duration,
      Simulated_Duration = T,
      Model = "Rapid"
    ))
    curves_ram[[length(curves_ram) + 1]] <- P
  }
}


# -------------------------------
# Gradual Abandonment Model 
# -------------------------------

for (i in 1:n_each_model) {
  T <- sample(duration_range[1]:duration_range[2], 1)
  t <- 0:(T - 1)

  # Sample parameters
  P0 <- runif(1, initial_pop_range[1], initial_pop_range[2])
  K <- runif(1, max(P0 + 1, carrying_capacity_range[1]), carrying_capacity_range[2])
  r1 <- runif(1, growth_rate_range_pos[1], growth_rate_range_pos[2])
  r2 <- runif(1, growth_rate_range_neg[1], growth_rate_range_neg[2])
  S <- runif(1, household_size_range[1], household_size_range[2])
  L <- rnorm(1, mean = mean(house_lifetime_range), sd = 5)
  if (L < house_lifetime_range[1] || L > house_lifetime_range[2]) next

  # Define durations for growth and decline
  grow_T <- round(rnorm(1, mean = T / 2, sd = 20))
  grow_T <- min(max(1, grow_T), T - 1)
  decline_T <- T - grow_T

  if (S <= 0 || L <= 0 || grow_T <= 0 || decline_T <= 0) next

  # Simulate two-phase population trajectory
  P_growth <- logistic_growth(P0, K, r1, grow_T)
  P_peak <- P_growth[length(P_growth)]
  P_decline <- P_peak * exp(r2 * (0:(decline_T - 1)))
  P <- c(P_growth, P_decline)
  max_pop <- max(P)

  if (abs(K - max_pop) >= 10) next

  # Estimate household counts and duration
  steps <- floor(T / L)
  house_counts <- sapply(0:(steps - 1), function(i) {
    index <- i * floor(L) + 1
    if (index <= length(P)) return(P[index] / S) else return(0)
  })
  H <- sum(unlist(house_counts))
  duration <- estimate_duration(P)


  # Enforce constraints: initial houses should not exceed 50 and last year should not have less than 8 inhabitants
  if (P[1] / S > 50 || P[length(P)] <= 0 || P[length(P)] >= 8) next
  
  # Final constraint: match observed archaeological household range
  if (H >= observed_min && H <= observed_max && duration <= duration_max) {
  # Write results   
	results_gam <- rbind(results_gam, data.frame(
      Initial_Population = P[1],
      Final_Population = P[length(P)],
      Max_Population = max_pop,
      Growth_Rate = r1,
      Decline_Rate = r2,
      Duration = duration,
      Simulated_Duration = T,
      Model = "Gradual"
    ))
    curves_gam[[length(curves_gam) + 1]] <- P
  }
}

# Combine both model outputs
results_all <- bind_rows(results_ram, results_gam)


# =============================================================================
# PLOTTING & EXPORTING STATISTICS
# =============================================================================


# -------------------------------
# Define output path for saving results
# -------------------------------
output_dir <- "FILEPATH"

# Generic plot-saving function
save_plot <- function(plot, name, width = 8, height = 6) {
  ggsave(file.path(output_dir, paste0(name, ".png")), plot = plot, width = width, height = height, dpi = 300)
  ggsave(file.path(output_dir, paste0(name, ".svg")), plot = plot, width = width, height = height)
}


# -------------------------------
# Create and save histograms for RAM
# -------------------------------
plot_rapid_init <- ggplot(results_ram, aes(x = Initial_Population)) +
  geom_histogram(bins = 40, fill = "grey70", color = "black") +
  labs(title = "Initial Population (Rapid Model)", x = "Initial Population", y = "Frequency") +
  theme_bw()

plot_rapid_growth <- ggplot(results_ram, aes(x = Growth_Rate)) +
  geom_histogram(bins = 40, fill = "grey70", color = "black") +
  labs(title = "Growth Rate (Rapid Model)", x = "Growth Rate", y = "Frequency") +
  theme_bw()

plot_rapid_max <- ggplot(results_ram, aes(x = Max_Population)) +
  geom_histogram(bins = 40, fill = "grey70", color = "black") +
  labs(title = "Max Population (Rapid Model)", x = "Max Population", y = "Frequency") +
  theme_bw()

plot_rapid_duration <- ggplot(results_ram, aes(x = Simulated_Duration)) +
  geom_histogram(bins = 40, fill = "grey70", color = "black") +
  labs(title = "Duration (Rapid Model)", x = "Years", y = "Frequency") +
  theme_bw()

save_plot(plot_rapid_init, "hist_initial_population_rapid")
save_plot(plot_rapid_growth, "hist_growth_rate_rapid")
save_plot(plot_rapid_max, "hist_max_population_rapid")
save_plot(plot_rapid_duration, "hist_duration_rapid")


# -------------------------------
# Create and save histograms for GAM
# -------------------------------
plot_gradual_init <- ggplot(results_gam, aes(x = Initial_Population)) +
  geom_histogram(bins = 40, fill = "grey70", color = "black") +
  labs(title = "Initial Population (Gradual Model)", x = "Initial Population", y = "Frequency") +
  theme_bw()

plot_gradual_growth <- ggplot(results_gam, aes(x = Growth_Rate)) +
  geom_histogram(bins = 40, fill = "grey70", color = "black") +
  labs(title = "Growth Rate (Gradual Model)", x = "Growth Rate", y = "Frequency") +
  theme_bw()

plot_gradual_decline <- ggplot(results_gam, aes(x = Decline_Rate)) +
  geom_histogram(bins = 40, fill = "grey70", color = "black") +
  labs(title = "Decline Rate (Gradual Model)", x = "Decline Rate", y = "Frequency") +
  theme_bw()

plot_gradual_max <- ggplot(results_gam, aes(x = Max_Population)) +
  geom_histogram(bins = 40, fill = "grey70", color = "black") +
  labs(title = "Max Population (Gradual Model)", x = "Max Population", y = "Frequency") +
  theme_bw()

plot_gradual_duration <- ggplot(results_gam, aes(x = Simulated_Duration)) +
  geom_histogram(bins = 30, fill = "grey70", color = "black") +
  labs(title = "Duration (Gradual Model)", x = "Years", y = "Frequency") +
  theme_bw()

save_plot(plot_gradual_init, "hist_initial_population_gradual")
save_plot(plot_gradual_growth, "hist_growth_rate_gradual")
save_plot(plot_gradual_decline, "hist_decline_rate_gradual")
save_plot(plot_gradual_max, "hist_max_population_gradual")
save_plot(plot_gradual_duration, "hist_duration_gradual")


# -------------------------------
# Export Line Plots of Population Trajectories
# -------------------------------
library(tibble)

df_rapid <- do.call(rbind, lapply(seq_along(curves_ram[1:min(50, length(curves_ram))]), function(i) {
  tibble(Year = 1:length(curves_ram[[i]]),
         Population = curves_ram[[i]],
         Simulation = paste0("Sim_", i))
}))

plot_line_rapid <- ggplot(df_rapid, aes(x = Year, y = Population, group = Simulation)) +
  geom_line(color = "black", alpha = 0.6) +
  labs(title = "Population Curves: Rapid Model", x = "Year", y = "Population") +
  theme_bw()

df_gradual <- do.call(rbind, lapply(seq_along(curves_gam[1:min(50, length(curves_gam))]), function(i) {
  tibble(Year = 1:length(curves_gam[[i]]),
         Population = curves_gam[[i]],
         Simulation = paste0("Sim_", i))
}))

plot_line_gradual <- ggplot(df_gradual, aes(x = Year, y = Population, group = Simulation)) +
  geom_line(color = "black", alpha = 0.6) +
  labs(title = "Population Curves: Gradual Model", x = "Year", y = "Population") +
  theme_bw()

save_plot(plot_line_rapid, "line_population_curves_rapid")
save_plot(plot_line_gradual, "line_population_curves_gradual")


# -------------------------------
# Export 5-Number Summaries to Excel
# -------------------------------
summary_stats <- function(x) {
  fivenum_x <- fivenum(x)
  names(fivenum_x) <- c("Min", "Q1", "Median", "Q3", "Max")
  return(c(fivenum_x, Count = length(x)))
}

metrics <- c("Initial_Population", "Growth_Rate", "Decline_Rate", "Max_Population", "Simulated_Duration")
summary_df <- data.frame()

for (m in metrics) {
  for (mod in c("Rapid", "Gradual")) {
    values <- results_all %>% filter(Model == mod) %>% pull(!!sym(m))
    stat_row <- summary_stats(values)
    summary_df <- rbind(summary_df, data.frame(
      Metric = m,
      Model = mod,
      t(stat_row)
    ))
  }
}

output_path <- file.path(output_dir, "statistics_summary.xlsx")
write_xlsx(summary_df, output_path)
cat("Excel summary saved to:\n", output_path, "\n")