#!/usr/bin/env Rscript

# Reproducible analysis for "When Will Tibo Press the Button Again?"
# Uses only base R. Run from this directory: Rscript analysis.R

options(stringsAsFactors = FALSE)
tz <- "Australia/Sydney"
as_of <- as.POSIXct("2026-08-24 13:20:00", tz = tz)
window_days <- 50
window_start <- as_of - window_days * 86400

d <- read.csv("resets.csv", check.names = FALSE)
d$time <- as.POSIXct(d$timestamp_sydney, format = "%Y-%m-%d %H:%M:%S", tz = tz)
d <- d[order(d$time), ]
d$gap_days <- c(NA, as.numeric(diff(d$time), units = "days"))

# A deliberately modest model: resets are events from a homogeneous Poisson
# process over the trailing 50 days. Jeffreys' prior for the rate gives
# lambda | data ~ Gamma(k + 1/2, exposure). Integrating over lambda produces
# a Lomax posterior-predictive waiting-time distribution.
k <- sum(d$time >= window_start & d$time <= as_of)
alpha <- k + 0.5
beta <- window_days
qwait <- function(p) beta * ((1 - p)^(-1 / alpha) - 1)
pwait <- function(days) 1 - (beta / (beta + days))^alpha
q50 <- qwait(0.50)
q75 <- qwait(0.75)
q90 <- qwait(0.90)
rate_mean <- alpha / beta
rate_lo <- qgamma(0.10, shape = alpha, rate = beta)
rate_hi <- qgamma(0.90, shape = alpha, rate = beta)

launch <- as.POSIXct("2026-07-10 00:00:00", tz = tz)
first <- min(d$time)
pre_n <- sum(d$time < launch)
post_n <- sum(d$time >= launch)
pre_days <- as.numeric(difftime(launch, first, units = "days"))
post_days <- as.numeric(difftime(as_of, launch, units = "days"))
pre_rate <- pre_n / pre_days
post_rate <- post_n / post_days
rate_ratio <- post_rate / pre_rate

fmt_time <- function(x) format(x, "%a %d %b, %H:%M %Z", tz = tz)
json <- sprintf(paste0(
  "{\n",
  "  \"as_of\": \"%s\",\n",
  "  \"events_total\": %d,\n",
  "  \"events_in_trailing_window\": %d,\n",
  "  \"window_days\": %d,\n",
  "  \"posterior_mean_resets_per_day\": %.4f,\n",
  "  \"posterior_rate_80pct_interval\": [%.4f, %.4f],\n",
  "  \"median_wait_days\": %.3f,\n",
  "  \"median_next_reset\": \"%s\",\n",
  "  \"probability_within_3_days\": %.4f,\n",
  "  \"probability_within_7_days\": %.4f,\n",
  "  \"probability_within_14_days\": %.4f,\n",
  "  \"post_launch_rate_ratio\": %.3f\n",
  "}\n"),
  format(as_of, "%Y-%m-%dT%H:%M:%S%z", tz = tz), nrow(d), k, window_days,
  rate_mean, rate_lo, rate_hi, q50,
  fmt_time(as_of + q50 * 86400), pwait(3), pwait(7), pwait(14), rate_ratio)
writeLines(json, "forecast.json")

bg <- "#0d0c0a"; surface <- "#171512"; body <- "#c9c2b7"
muted <- "#827b70"; gold <- "#e3b562"; cyan <- "#63d8d2"; red <- "#ff746c"
chart_par <- function(mar = c(4.2, 4.6, 2.2, 1.0)) {
  par(bg = bg, fg = body, col.axis = muted, col.lab = body, col.main = "#f4efe7",
      family = "sans", mar = mar, las = 1, bty = "n", xaxs = "i", yaxs = "i")
}
grid_y <- function(vals) abline(h = vals, col = "#29251f", lwd = 1)

# 1. Cumulative resets
svg("cumulative-resets.svg", width = 11, height = 5.2, bg = bg, pointsize = 12)
chart_par(c(4.5, 4.3, 2.6, 1.2))
x <- as.Date(d$time, tz = tz)
plot_y_max <- ceiling((nrow(d) + 2) / 5) * 5
plot(range(x) + c(-2, 8), c(0, plot_y_max), type = "n", axes = FALSE,
     xlab = "2026", ylab = "Confirmed broad resets", main = "The button is getting warmer")
grid_y(seq(0, plot_y_max, 5))
axis(2, at = seq(0, plot_y_max, 5), col = NA, col.ticks = muted)
axis.Date(1, at = seq(as.Date("2026-05-17"), as.Date(as_of, tz = tz), by = "2 weeks"),
          format = "%d %b", col = NA, col.ticks = muted)
lines(x, seq_len(nrow(d)), type = "s", lwd = 3, col = gold)
points(x, seq_len(nrow(d)), pch = 21, cex = 1.15, bg = bg, col = gold, lwd = 2)
abline(v = as.Date("2026-07-10"), col = cyan, lty = 2, lwd = 1.5)
text(as.Date("2026-07-10") + 1, 3.1, "ChatGPT Work / Sol launch", col = cyan, adj = 0, cex = .85)
text(max(x), nrow(d) + .7, as.character(nrow(d)), col = gold, font = 2, cex = 1.25)
dev.off()

# 2. Gaps
svg("reset-gaps.svg", width = 11, height = 5.2, bg = bg, pointsize = 12)
chart_par(c(5.2, 4.3, 2.6, 1.2))
g <- d$gap_days[-1]
labels <- format(as.Date(d$time[-1], tz = tz), "%d %b")
cols <- ifelse(as.Date(d$time[-1], tz = tz) >= as.Date("2026-07-10"), gold, muted)
plot(seq_along(g), g, type = "n", ylim = c(0, max(g) * 1.12), axes = FALSE,
     xlab = "Reset date", ylab = "Days since previous reset", main = "Time between resets")
grid_y(seq(0, 14, 2))
segments(seq_along(g), 0, seq_along(g), g, col = paste0(cols, "99"), lwd = 4)
points(seq_along(g), g, pch = 21, bg = bg, col = cols, lwd = 2, cex = 1.15)
axis(2, at = seq(0, 14, 2), col = NA, col.ticks = muted)
axis(1, at = seq_along(g), labels = labels, las = 2, cex.axis = .72, col = NA, col.ticks = muted)
abline(h = median(g), col = cyan, lty = 2, lwd = 1.5)
text(1, median(g) + .55, sprintf("median gap: %.1f days", median(g)), col = cyan, adj = 0, cex = .85)
legend("topright", legend = c("Before 10 Jul", "10 Jul onward"), col = c(muted, gold),
       pch = 16, bty = "n", text.col = body, cex = .85)
dev.off()

# 3. Posterior-predictive forecast
svg("next-reset-forecast.svg", width = 11, height = 5.5, bg = bg, pointsize = 12)
chart_par(c(4.8, 4.6, 2.8, 1.2))
t <- seq(0, 14, length.out = 400)
p <- pwait(t)
plot(t, p * 100, type = "n", ylim = c(0, 100), axes = FALSE,
     xlab = sprintf("Days after %s", format(as.Date(as_of, tz = tz), "%d %B")), ylab = "Chance at least one reset has happened",
     main = "The next press: posterior-predictive probability")
grid_y(seq(0, 100, 20))
polygon(c(t, rev(t)), c(rep(0, length(t)), rev(p * 100)), col = "#e3b56222", border = NA)
lines(t, p * 100, col = gold, lwd = 4)
axis(2, at = seq(0, 100, 20), labels = paste0(seq(0, 100, 20), "%"), col = NA, col.ticks = muted)
axis(1, at = c(0, 2, 4, 6, 8, 10, 12, 14),
     labels = format(as.Date(as_of, tz = tz) + c(0, 2, 4, 6, 8, 10, 12, 14), "%d %b"),
     col = NA, col.ticks = muted)
abline(h = 50, col = "#ffffff33", lty = 3)
abline(v = q50, col = cyan, lty = 2, lwd = 2)
points(c(3, 7, 14), pwait(c(3, 7, 14)) * 100, pch = 21, bg = bg, col = gold, lwd = 2, cex = 1.25)
text(q50 + .25, 44, sprintf("median: %s", format(as.Date(as_of + q50 * 86400), "%d %b")),
     col = cyan, adj = 0, cex = .9, font = 2)
for (z in c(3, 7, 14)) text(z, pwait(z) * 100 + 7,
  sprintf("%.0f%%", 100 * pwait(z)), col = gold, cex = .9, font = 2)
dev.off()

# 4. Social card (1200 x 630)
png("social-card.png", width = 1200, height = 630, bg = bg, res = 144)
par(bg = bg, fg = body, family = "sans", mar = rep(0, 4), xaxs = "i", yaxs = "i")
plot.new(); plot.window(xlim = c(0, 1200), ylim = c(0, 630))
rect(0, 0, 1200, 630, col = bg, border = NA)
rect(0, 0, 18, 630, col = gold, border = NA)
text(72, 550, "THE TIBO RESET INDEX", adj = 0, col = cyan, cex = 1.15, font = 2)
text(72, 455, "When will he press", adj = 0, col = "#f4efe7", cex = 2.7, font = 2)
text(72, 372, "the button again?", adj = 0, col = "#f4efe7", cex = 2.7, font = 2)
observation_days <- as.integer(as.Date(as_of, tz = tz) - as.Date(first, tz = tz)) + 1
social_horizon <- 3
social_probability <- round(100 * pwait(social_horizon))
text(72, 260, as.character(nrow(d)), adj = 0, col = gold, cex = 4.4, font = 2)
text(236, 267, sprintf("broad resets\nin %d days", observation_days), adj = 0, col = body, cex = 1.35, font = 2)
text(620, 260, sprintf("%d%%", social_probability), adj = 0, col = gold, cex = 4.4, font = 2)
text(888, 267, sprintf("chance of another\nwithin %d days", social_horizon), adj = 0, col = body, cex = 1.35, font = 2)
segments(72, 148, 1128, 148, col = "#312d27", lwd = 2)
text(72, 88, sprintf("R forecast: median next reset · %s", fmt_time(as_of + q50 * 86400)), adj = 0,
     col = "#d9d2c8", cex = 1.13)
text(1128, 35, "wasnotwas.com", adj = 1, col = muted, cex = .9)
dev.off()

cat(sprintf("%d resets total; %d in trailing %d days\n", nrow(d), k, window_days))
cat(sprintf("Median next reset: %s (%.2f days)\n", fmt_time(as_of + q50 * 86400), q50))
cat(sprintf("P(within 3/7/14 days): %.1f%% / %.1f%% / %.1f%%\n",
            100*pwait(3), 100*pwait(7), 100*pwait(14)))
cat(sprintf("Post-launch descriptive rate ratio: %.2fx\n", rate_ratio))