# =============================================================================
# generate_agentic_reorder_images.R
# "I Gave an AI Agent the Reorder Button" — June 2026 flagship post
# -----------------------------------------------------------------------------
# Thesis under test: the danger of an agentic AI replenishment system is not
# that it is dumb — it is that it is TOO RESPONSIVE. By chasing recent demand
# it manufactures its own bullwhip, amplifying demand variance into order
# variance. The fix is guardrails (forecast smoothing, dampened correction,
# order bounds), not more autonomy.
#
# We simulate ONE representative SKU over 250 days under three replenishment
# policies that drive the reorder decision, then measure the bullwhip ratio,
# service level, inventory, and cost for each. Numbers emerge from the
# simulation — nothing is hand-tuned to flatter the thesis.
#
# Run from project root:  Rscript Scripts/generate_agentic_reorder_images.R
# =============================================================================

source("Scripts/theme_inphronesys.R")

suppressPackageStartupMessages({
  library(ggplot2)
  library(dplyr)
  library(tidyr)
  library(scales)
  library(patchwork)
  library(jsonlite)
  library(ggrepel)
})

set.seed(42)

img_dir  <- "Images"
dash_dir <- "Dashboards"
if (!dir.exists(img_dir))  dir.create(img_dir)
if (!dir.exists(dash_dir)) dir.create(dash_dir)

# =============================================================================
# 1. SCENARIO PARAMETERS  (the single source of truth)
# =============================================================================
N_DAYS        <- 250      # simulation horizon (days)
WARMUP        <- 30       # days excluded from metrics (transient settle)
LEAD_TIME     <- 5        # fixed replenishment lead time (days)
BASE_DEMAND   <- 100      # baseline mean daily demand (units)
SEASON_AMP    <- 12       # weekly seasonality amplitude (units)
NOISE_SD      <- 14       # daily Gaussian demand noise (units)
STEP_DAY      <- 150      # day a permanent demand regime shift hits
STEP_SIZE     <- 30       # size of the level shift (units/day)

# Costs (per unit unless noted) — lost-sales model (NOT backorder)
HOLD_COST     <- 0.50     # holding cost per unit per day on hand
ORDER_COST    <- 40       # fixed cost per replenishment order placed
STOCKOUT_PEN  <- 8        # lost-margin penalty per unit of demand lost
Z_SERVICE     <- 1.65     # safety factor (~95% cycle service)

EVAL <- (WARMUP + 1):N_DAYS    # evaluation window for all metrics

# =============================================================================
# 2. DEMAND GENERATION
# =============================================================================
# Daily demand = base level + permanent step + weekly seasonality + noise,
# floored at zero. Lost sales: unmet demand is lost, never backordered.
gen_demand <- function() {
  t      <- 1:N_DAYS
  level  <- BASE_DEMAND + ifelse(t >= STEP_DAY, STEP_SIZE, 0)
  season <- SEASON_AMP * sin(2 * pi * t / 7)
  noise  <- rnorm(N_DAYS, 0, NOISE_SD)
  pmax(0, round(level + season + noise))
}
demand <- gen_demand()

# =============================================================================
# 3. INVENTORY SIMULATION ENGINE  (shared across all three policies)
# =============================================================================
# Daily sequence:
#   (1) receive the order that was placed LEAD_TIME days ago
#   (2) demand arrives; sales = min(demand, on-hand); shortfall is LOST
#   (3) inventory position IP = on-hand + on-order (in transit)
#   (4) the policy decides today's order, which arrives in LEAD_TIME days
simulate_policy <- function(demand, policy_fn, L = LEAD_TIME) {
  n        <- length(demand)
  on_hand  <- numeric(n)
  orders   <- numeric(n)
  sales    <- numeric(n)
  lost     <- numeric(n)
  ip_rec   <- numeric(n)
  pipeline <- numeric(n + L + 2)   # arrivals indexed by arrival day

  # Steady-state warm start: stock on hand at order-up-to level and seed the
  # pipeline so the first L days have replenishment already in transit.
  S0  <- BASE_DEMAND * (L + 1) + Z_SERVICE * NOISE_SD * sqrt(L + 1)
  oh  <- S0
  for (d in 1:L) pipeline[d] <- BASE_DEMAND

  state <- list()
  for (t in 1:n) {
    oh <- oh + pipeline[t]                 # (1) receive
    d  <- demand[t]                        # (2) demand & lost sales
    s  <- min(d, oh)
    oh <- oh - s
    l  <- d - s
    on_order <- sum(pipeline[(t + 1):(t + L)])
    IP <- oh + on_order                    # (3) inventory position

    res   <- policy_fn(t, demand, IP, state, L)   # (4) decide order
    o     <- max(0, res$order)
    state <- res$state
    pipeline[t + L] <- pipeline[t + L] + o

    on_hand[t] <- oh
    orders[t]  <- o
    sales[t]   <- s
    lost[t]    <- l
    ip_rec[t]  <- IP
  }

  data.frame(day = 1:n, demand = demand, order = orders,
             on_hand = on_hand, sales = sales, lost = lost, ip = ip_rec)
}

# =============================================================================
# 4. THE THREE REPLENISHMENT POLICIES
# =============================================================================

# --- Policy 1: Static reorder-point planner (the human baseline) -------------
# Order-up-to (base-stock) with a FIXED mean and sigma set once. Does not
# react to daily swings. Orders therefore track demand 1:1 (bullwhip ~ 1),
# but the static target is blind to the regime shift at day 150.
S_STATIC <- BASE_DEMAND * (LEAD_TIME + 1) + Z_SERVICE * NOISE_SD * sqrt(LEAD_TIME + 1)
policy_static <- function(t, demand, IP, state, L) {
  list(order = S_STATIC - IP, state = state)
}

# --- Policy 2: Naive reactive agent (the "give the LLM the button" case) -----
# Base-stock target driven by a short 3-day moving-average forecast and a
# short-window sigma. No dampening, no bounds. The (L+1) multiplier on a
# jumpy forecast is exactly the Chen et al. (2000) amplification mechanism.
NAIVE_P   <- 3      # forecast window (days) — very reactive
NAIVE_SIG <- 10     # sigma estimation window (days)
policy_naive <- function(t, demand, IP, state, L) {
  hist  <- demand[max(1, t - NAIVE_P):(t - 1)]
  if (t == 1) hist <- demand[1]
  Dhat  <- mean(hist)
  sig   <- if (t > NAIVE_SIG) sd(demand[(t - NAIVE_SIG):(t - 1)]) else NOISE_SD
  S     <- Dhat * (L + 1) + Z_SERVICE * sig * sqrt(L + 1)
  list(order = S - IP, state = state)
}

# --- Policy 3: Guardrailed agent (forecast-smoothed, dampened, bounded) ------
# It chases the SAME base-stock target as the naive agent, but with three
# guardrails: (a) an exponentially-smoothed demand forecast (alpha = 0.15)
# so the target itself does not jump; (b) a dampened order signal — the order
# only moves a fraction (beta = 0.25) of the way toward the raw correction
# each day, the control-theoretic proportional dampening of Dejonckheere &
# Disney (2003); (c) a hard order bound at 1.8x the smoothed daily forecast.
# Same autonomy, same target, far less amplification.
GUARD_ALPHA <- 0.15   # forecast smoothing constant
GUARD_BETA  <- 0.25   # order dampening (fraction of correction applied per day)
GUARD_SIGW  <- 21     # sigma estimation window (days, ~3 weeks)
GUARD_CAP   <- 1.8    # order ceiling as multiple of smoothed daily forecast
policy_guard <- function(t, demand, IP, state, L) {
  if (is.null(state$dhat)) { state$dhat <- demand[1]; state$last_order <- demand[1] }
  if (t > 1) state$dhat <- GUARD_ALPHA * demand[t - 1] + (1 - GUARD_ALPHA) * state$dhat
  Dhat   <- state$dhat
  sig    <- if (t > GUARD_SIGW) sd(demand[(t - GUARD_SIGW):(t - 1)]) else NOISE_SD
  S      <- Dhat * (L + 1) + Z_SERVICE * sig * sqrt(L + 1)
  target <- S - IP                                       # the raw base-stock order
  raw    <- (1 - GUARD_BETA) * state$last_order + GUARD_BETA * target  # dampened
  ord    <- min(max(0, raw), GUARD_CAP * Dhat)           # order bound
  state$last_order <- ord
  list(order = ord, state = state)
}

# =============================================================================
# 5. RUN THE SIMULATIONS
# =============================================================================
sim_static <- simulate_policy(demand, policy_static)
sim_naive  <- simulate_policy(demand, policy_naive)
sim_guard  <- simulate_policy(demand, policy_guard)

# =============================================================================
# 6. METRICS
# =============================================================================
bullwhip_ratio <- function(orders, dmd) var(orders) / var(dmd)

compute_metrics <- function(sim, label) {
  e        <- sim[EVAL, ]
  fill     <- sum(e$sales) / sum(e$demand)             # units fill rate
  avg_inv  <- mean(e$on_hand)
  so_days  <- sum(e$lost > 0)
  lost_u   <- sum(e$lost)
  n_orders <- sum(e$order > 0.5)
  hold_c   <- HOLD_COST * sum(e$on_hand)
  order_c  <- ORDER_COST * n_orders
  pen_c    <- STOCKOUT_PEN * lost_u
  total_c  <- hold_c + order_c + pen_c
  bw       <- bullwhip_ratio(e$order, e$demand)
  data.frame(
    policy        = label,
    fill_rate     = fill,
    avg_inventory = avg_inv,
    stockout_days = so_days,
    lost_units    = lost_u,
    n_orders      = n_orders,
    holding_cost  = hold_c,
    ordering_cost = order_c,
    stockout_cost = pen_c,
    total_cost    = total_c,
    bullwhip      = bw,
    stringsAsFactors = FALSE
  )
}

metrics <- bind_rows(
  compute_metrics(sim_static, "Static planner"),
  compute_metrics(sim_naive,  "Naive agent"),
  compute_metrics(sim_guard,  "Guardrailed agent")
)

policy_levels <- c("Static planner", "Naive agent", "Guardrailed agent")
policy_cols   <- c("Static planner"    = iph_colors$grey,
                   "Naive agent"       = iph_colors$red,
                   "Guardrailed agent" = iph_colors$blue)

# --- Chart 1: orders vs demand across the three policies (amplification) ------
ts_df <- bind_rows(
  transform(sim_static, policy = "Static planner"),
  transform(sim_naive,  policy = "Naive agent"),
  transform(sim_guard,  policy = "Guardrailed agent")
) %>%
  filter(day >= WARMUP + 1) %>%
  mutate(policy = factor(policy, levels = policy_levels))

ts_long <- ts_df %>%
  select(day, policy, Demand = demand, Orders = order) %>%
  pivot_longer(c(Demand, Orders), names_to = "series", values_to = "units")

p1 <- ggplot(ts_long, aes(x = day, y = units, color = series, linewidth = series)) +
  geom_vline(xintercept = STEP_DAY, linetype = "dashed",
             color = iph_colors$dark, alpha = 0.4) +
  geom_line() +
  facet_wrap(~ policy, ncol = 1) +
  scale_color_manual(values = c("Demand" = iph_colors$lightgrey,
                                "Orders" = iph_colors$navy)) +
  scale_linewidth_manual(values = c("Demand" = 0.9, "Orders" = 0.6), guide = "none") +
  labs(title = "Orders chase demand, and the naive agent overshoots",
       subtitle = "Daily orders (navy) vs demand (grey). Dashed line = the +30 unit/day demand shift",
       x = "Day", y = "Units", color = NULL) +
  theme_inphronesys(grid = "y") +
  theme(legend.position = "bottom")

ggsave(file.path(img_dir, "agent_orders_vs_demand.png"), p1,
       width = 8, height = 7, dpi = 100, bg = "white")

# --- Chart 2: bullwhip ratio bar chart (HERO) --------------------------------
bw_df <- metrics %>% mutate(policy = factor(policy, levels = policy_levels))

p2 <- ggplot(bw_df, aes(x = policy, y = bullwhip, fill = policy)) +
  geom_col(width = 0.62, show.legend = FALSE) +
  geom_hline(yintercept = 1, linetype = "dashed", color = iph_colors$dark) +
  geom_text(aes(label = sprintf("%.2fx", bullwhip)),
            vjust = -0.5, fontface = "bold", size = 4.6, color = iph_colors$dark) +
  scale_fill_manual(values = policy_cols) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.16))) +
  labs(title = "The bullwhip ratio: who amplifies demand variance?",
       subtitle = "Var(orders) / Var(demand). Above 1.0, the policy amplifies the demand it serves",
       x = NULL, y = "Bullwhip ratio") +
  theme_inphronesys(grid = "y")

ggsave(file.path(img_dir, "agent_bullwhip_ratio.png"), p2,
       width = 8, height = 5, dpi = 100, bg = "white")

# --- Chart 3: service level vs average inventory trade-off -------------------
trade_df <- metrics %>%
  mutate(policy   = factor(policy, levels = policy_levels),
         fill_pct = fill_rate * 100,
         lab      = sprintf("%s\nbullwhip %.2fx", policy, bullwhip))

p3 <- ggplot(trade_df, aes(x = avg_inventory, y = fill_pct)) +
  geom_point(aes(size = bullwhip, color = policy), alpha = 0.78) +
  geom_text_repel(aes(label = lab, color = policy),
                  fontface = "bold", size = 3.6, lineheight = 0.95,
                  segment.color = iph_colors$grey, show.legend = FALSE) +
  scale_color_manual(values = policy_cols, guide = "none") +
  scale_size_continuous(range = c(5, 19), name = "Bullwhip ratio") +
  labs(title = "Same service, same stock: the only difference is the bullwhip",
       subtitle = "Up-and-left is better: higher fill, less inventory. Bubble size = bullwhip ratio",
       x = "Average on-hand inventory (units)", y = "Fill rate") +
  theme_inphronesys(grid = "xy") +
  theme(legend.position = "bottom")

ggsave(file.path(img_dir, "agent_service_inventory_tradeoff.png"), p3,
       width = 8, height = 5, dpi = 100, bg = "white")