# MSC-P-048 candidate (local, unpublished): reading an A/B test before its planned date.
# Base R only. Same declared recursion, design, random per-visitor assignment, binomial draw,
# decision rules, boundary, invariant checks and exact size as msc-p048-reference.py; the two
# programs must print the same lines. Shares are printed as a count and five decimals: a rate out
# of 20000 is exact at five decimals, while four decimals land on half-way values that the two
# languages round in opposite directions.
LCG_SEED <- 20260916
LCG_MODULUS <- 2147483648
LCG_INCREMENT <- 12345
LCG_HIGH <- 16838   # 1103515245 = 16838 * 65536 + 20077, split to keep products exact
LCG_LOW <- 20077

BASE <- 0.05
UPLIFT <- 0.015
PER_DAY <- 250L  # visitors a day in total, each assigned at random to one arm
PLANNED_DAYS <- 28L
MONITOR_DAYS <- 84L
ALPHA <- 0.05
Z_FIXED <- 1.959963985
LOOK_POINTS <- c(7L, 14L, 28L, 84L)
RHO_DAYS <- c(7L, 28L, 56L)  # calibration dates declared a priori, beside the source optimum
REPLICATIONS <- 20000L
EXAMPLE_LIMIT <- 200L
GRID_STEPS <- 2000L  # grid on which the sub-Gaussian inequality is checked
PROBE_DAYS <- 2000L  # days drawn by the invariants to measure the assignment and the coin flip
WALKS <- 5000L
WALK_STEPS <- 2000L
GUARD_REPLICATIONS <- 2000L  # replications of the blocking null-hypothesis check
STICKY_REPLICATIONS <- 5000L  # replications of the sticky-assignment measurement
STICKY_BLOCK <- 10L  # visitors sharing one arm when the assignment is drawn but sticky

lcg_state <- LCG_SEED
uniform <- function() {
  high <- (LCG_HIGH * lcg_state) %% LCG_MODULUS
  lcg_state <<- (high * 65536 + LCG_LOW * lcg_state + LCG_INCREMENT) %% LCG_MODULUS
  (lcg_state + 0.5) / LCG_MODULUS
}

# (1 - p) to the n for n = 0 to limit, by repeated multiplication: the same operation in both
# languages, rather than a power, whose last bit differs between Python and R.
first_terms <- function(p, limit) {
  terms <- numeric(limit + 1L)
  terms[1] <- 1
  for (i in seq_len(limit)) terms[i + 1L] <- terms[i] * (1 - p)
  terms
}

TERMS_HALF <- first_terms(0.5, PER_DAY)
TERMS_BASE <- first_terms(BASE, PER_DAY)
TERMS_UPLIFT <- first_terms(BASE + UPLIFT, PER_DAY)

# Exact binomial draw by inversion: one uniform, then the probability recursion.
binomial <- function(n, p, terms) {
  u <- uniform()
  ratio <- p / (1 - p)
  term <- terms[n + 1L]
  total <- term
  k <- 0L
  while (u > total && k < n) {
    term <- term * (n - k) / (k + 1) * ratio
    k <- k + 1L
    total <- total + term
  }
  k
}

# Two-sided pooled z statistic of two proportions, arms of possibly different sizes.
pooled_z <- function(n_control, control, n_variant, variant) {
  if (n_control <= 0 || n_variant <= 0) return(0)
  pooled <- (control + variant) / (n_control + n_variant)
  if (pooled <= 0 || pooled >= 1) return(0)
  spread <- sqrt(pooled * (1 - pooled) * (1 / n_control + 1 / n_variant))
  (variant / n_variant - control / n_control) / spread
}

# Two-sided normal-mixture boundary on the conversion gap, read on the conversion clock. Closed
# form of the two-sided normal mixture, equation (14) of Howard, Ramdas, McAuliffe and Sekhon
# (2018), with the initial likelihood ratio l0 = 1. Proposition 6 of the same article is the
# one-sided form and is not used here.
mixture_bound <- function(conversions, rho) {
  shifted <- conversions + rho
  sqrt(shifted * log(shifted / (rho * ALPHA * ALPHA)))
}

# Rho that minimises the threshold at a target number of conversions, by ternary search. The
# source states the boundary is tightest at a fixed multiple of rho, its equation (21); the search
# reproduces that multiple instead of hard-coding it, and an invariant checks the ratio.
tightest_rho <- function(target) {
  low <- 0.1
  high <- target
  for (i in seq_len(200L)) {
    first <- low + (high - low) / 3
    second <- high - (high - low) / 3
    if (mixture_bound(target, first) < mixture_bound(target, second)) high <- second else low <- first
  }
  (low + high) / 2
}

# Conversions expected by a date, under the declared base rate and daily traffic.
conversions_expected <- function(day) PER_DAY * day * BASE

# A day where the assignment is drawn but sticky: blocks of visitors share one arm. This is what
# cookie-level or session-level assignment does in a real tool. The shares stay at one half on
# average, so the measured invariants pass, but the conversion labels are no longer independent
# and the boundary loses its guarantee. The programs measure how much.
draw_sticky_day <- function(uplift) {
  n_variant <- 0L
  blocks <- PER_DAY %/% STICKY_BLOCK
  for (b in seq_len(blocks)) if (uniform() < 0.5) n_variant <- n_variant + STICKY_BLOCK
  n_control <- PER_DAY - n_variant
  terms_variant <- if (uplift > 0) TERMS_UPLIFT else TERMS_BASE
  variant <- binomial(n_variant, BASE + uplift, terms_variant)
  control <- binomial(n_control, BASE, TERMS_BASE)
  c(n_control, control, n_variant, variant)
}

# One day: the split of visitors between arms, then the conversions of each arm.
draw_day <- function(uplift) {
  n_variant <- binomial(PER_DAY, 0.5, TERMS_HALF)
  n_control <- PER_DAY - n_variant
  terms_variant <- if (uplift > 0) TERMS_UPLIFT else TERMS_BASE
  variant <- binomial(n_variant, BASE + uplift, terms_variant)
  control <- binomial(n_control, BASE, TERMS_BASE)
  c(n_control, control, n_variant, variant)
}

# Checks run before anything is computed; any failure stops the run. Two of them are measured on
# the drawn data, the assignment and the conversion coin flip, and can fail. The sub-Gaussian
# inequality is a mathematical identity, verified numerically rather than assumed; its least
# favourable point is the first of the grid.
check_invariants <- function() {
  failures <- character(0)
  worst <- NA_real_
  for (i in seq_len(GRID_STEPS)) {
    lam <- i * 0.01
    gap <- log(cosh(lam)) - lam * lam / 2
    worst <- if (is.na(worst) || gap > worst) gap else worst
  }
  if (worst > 0) failures <- c(failures, "sub_gaussian_inequality")

  visitors <- 0
  to_variant <- 0
  conversions <- 0
  variant_conversions <- 0
  for (i in seq_len(PROBE_DAYS)) {
    day <- draw_day(0)
    visitors <- visitors + day[1] + day[3]
    to_variant <- to_variant + day[3]
    conversions <- conversions + day[2] + day[4]
    variant_conversions <- variant_conversions + day[4]
  }
  assignment <- to_variant / visitors
  coin <- variant_conversions / conversions
  if (abs(assignment - 0.5) > 0.005) failures <- c(failures, "balanced_assignment")
  if (abs(coin - 0.5) > 0.01) failures <- c(failures, "conversion_coin_flip")

  target <- conversions_expected(PLANNED_DAYS)
  rho <- tightest_rho(target)
  ratio <- if (rho > 0) target / rho else 0
  if (rho <= 0) failures <- c(failures, "positive_rho")
  if (abs(ratio - 8.212) > 0.01) failures <- c(failures, "tightest_point_matches_source")
  for (probe in c(1, 10, 100, 500, 2000, target, 10 * target)) {
    shifted <- probe + rho
    rearranged <- shifted * (log(shifted) - log(rho) - 2 * log(ALPHA))
    if (abs(mixture_bound(probe, rho)^2 - rearranged) > 1e-6) {
      failures <- c(failures, "boundary_closed_form")
      break
    }
  }

  saved <- lcg_state
  lcg_state <<- LCG_SEED
  before <- lcg_state
  invisible(binomial(PER_DAY, BASE, TERMS_BASE))
  # Parenthesised like the arm counts below: %% binds tighter than *, so make the order explicit.
  expected <- (((LCG_HIGH * before) %% LCG_MODULUS) * 65536 + LCG_LOW * before + LCG_INCREMENT) %% LCG_MODULUS
  if (lcg_state != expected) failures <- c(failures, "one_uniform_per_draw")

  # The clock the boundary is actually read on, checked end to end on a short test: the verdicts
  # must be reproducible from the conversion total alone, recomputed outside run_test.
  lcg_state <<- LCG_SEED
  looks <- run_test(30L, 0, rho)
  lcg_state <<- LCG_SEED
  control <- 0
  variant <- 0
  for (day in seq_len(30L)) {
    drawn <- draw_day(0)
    control <- control + drawn[2]
    variant <- variant + drawn[4]
    if (looks$gap[day] != variant - control) {
      failures <- c(failures, "boundary_reads_conversion_gap")
      break
    }
    clock_conversions <- control + variant
    expected_cross <- if (clock_conversions > 0) abs(variant - control) >= mixture_bound(clock_conversions, rho) else FALSE
    if (looks$crossed[day, 1] != expected_cross) {
      failures <- c(failures, "boundary_reads_conversion_clock")
      break
    }
  }
  lcg_state <<- saved

  cat(sprintf("metric.invariant.sub_gaussian_worst_gap=%.3e grid=%d_values_lambda_0.01_to_%.2f\n",
              worst, GRID_STEPS, GRID_STEPS * 0.01))
  cat(sprintf("metric.invariant.measured visitors=%d to_variant_share=%.5f conversions=%d variant_share=%.5f\n",
              visitors, assignment, conversions, coin))
  cat(sprintf("metric.invariant.tightest_point target_over_rho=%.3f source_states=8.212\n", ratio))
  cat(sprintf("metric.invariant.checks_failed=%d %s\n", length(failures),
              if (length(failures)) paste(failures, collapse = ",") else "none"))
  if (length(failures)) stop(paste0("invariant failed: ", paste(failures, collapse = ",")))
}

# Independent check of the boundary: a Rademacher walk, where the null holds by construction.
boundary_crossing_check <- function(walks, steps, rho) {
  crossed <- 0L
  for (w in seq_len(walks)) {
    position <- 0L
    for (step in seq_len(steps)) {
      position <- position + if (uniform() < 0.5) 1L else -1L
      if (abs(position) >= mixture_bound(step, rho)) {
        crossed <- crossed + 1L
        break
      }
    }
  }
  crossed
}

log_choose <- function(n, k) lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1)

# Exact two-sided size of the pooled z test with n per arm, with no simulation. For a fixed control
# count c the test rejects when 2n (v - c)^2 >= z^2 (c + v)(2n - c - v), a quadratic in v whose
# outside is the rejection set; the exact probability then follows from the distribution function.
exact_size <- function(n) {
  k_all <- 0:n
  pmf <- exp(log_choose(n, k_all) + k_all * log(BASE) + (n - k_all) * log(1 - BASE))
  cdf <- cumsum(pmf)
  k <- Z_FIXED * Z_FIXED
  size <- 0
  a <- 2 * n + k
  for (c in 0:n) {
    b <- -2 * (2 * n * c + k * n - k * c)
    constant <- 2 * n * c * c - k * (2 * n * c - c * c)
    disc <- b * b - 4 * a * constant
    if (disc <= 0) next
    root <- sqrt(disc)
    below <- floor((-b - root) / (2 * a))
    above <- ceiling((-b + root) / (2 * a))
    mass <- 0
    if (below >= 0) {
      mass <- mass + cdf[min(below, n) + 1]
      if (c == 0) mass <- mass - pmf[1]
    }
    if (above <= n) {
      index <- max(above, 0)
      mass <- mass + 1 - (if (index > 0) cdf[index] else 0)
    }
    size <- size + pmf[c + 1] * mass
  }
  size
}

# One experiment read every day: z statistic, conversion gap, and each boundary's verdict.
run_test <- function(days, uplift, rhos) {
  control <- 0L
  variant <- 0L
  n_control <- 0L
  n_variant <- 0L
  z <- numeric(days)
  gap <- integer(days)
  crossed <- matrix(FALSE, days, length(rhos))
  for (day in seq_len(days)) {
    drawn <- draw_day(uplift)
    n_control <- n_control + drawn[1]
    control <- control + drawn[2]
    n_variant <- n_variant + drawn[3]
    variant <- variant + drawn[4]
    z[day] <- pooled_z(n_control, control, n_variant, variant)
    gap[day] <- variant - control
    conversions <- control + variant
    if (conversions > 0) crossed[day, ] <- abs(variant - control) >= mixture_bound(conversions, rhos)
  }
  list(z = z, gap = gap, crossed = crossed)
}

summary_share <- function(hits, total) {
  share <- hits / total
  c(hits, share, sqrt(max(share * (1 - share), 0) / total))
}

first_day <- function(flags, limit) {
  hits <- which(flags[seq_len(limit)])
  if (length(hits)) hits[1] else 0L
}

check_invariants()

target <- conversions_expected(PLANNED_DAYS)
labels <- c("optimal", sprintf("day_%d", RHO_DAYS))
rhos <- c(tightest_rho(target), conversions_expected(RHO_DAYS))
reference <- 1L  # the source-optimal tuning is this page's reference

cat(sprintf("metric.design.base_rate=%.3f uplift_synthetic=%.3f visitors_per_day=%d assignment=random_half\n",
            BASE, UPLIFT, PER_DAY))
cat(sprintf("metric.design.planned_days=%d monitored_days=%d alpha=%.3f replications=%d\n",
            PLANNED_DAYS, MONITOR_DAYS, ALPHA, REPLICATIONS))
cat(sprintf("metric.design.visitors_expected_per_arm_at_planned_day=%.1f conversions_expected=%.1f\n",
            PER_DAY * PLANNED_DAYS / 2, target))
cat(sprintf("metric.design.z_fixed=%.6f boundary_clock=conversions\n", Z_FIXED))
for (j in seq_along(labels)) {
  over_root <- mixture_bound(target, rhos[j]) / sqrt(target)
  cat(sprintf("metric.design.boundary tuning=%s rho=%.4f threshold_over_sqrt_conversions=%.4f threshold_in_z=%.4f\n",
              labels[j], rhos[j], over_root, over_root / sqrt(1 - BASE)))
}

saved_state <- lcg_state
lcg_state <<- LCG_SEED
crossed_walks <- boundary_crossing_check(WALKS, WALK_STEPS, rhos[reference])
lcg_state <<- saved_state
walk_rate <- crossed_walks / WALKS
cat(sprintf("metric.check.boundary_on_rademacher_walk walks=%d steps=%d crossed=%d rate=%.5f alpha=%.3f verdict=%s\n",
            WALKS, WALK_STEPS, crossed_walks, walk_rate, ALPHA, if (walk_rate <= ALPHA) "pass" else "fail"))
if (walk_rate > ALPHA) stop("boundary check failed: crossing rate above alpha")

lcg_state <<- LCG_SEED
guard <- 0L
for (iter in seq_len(GUARD_REPLICATIONS)) {
  looks <- run_test(MONITOR_DAYS, 0, rhos[reference])
  if (first_day(looks$crossed[, 1], MONITOR_DAYS) > 0L) guard <- guard + 1L
}
guard_rate <- guard / GUARD_REPLICATIONS
cat(sprintf("metric.check.null_false_alarm replications=%d looks=%d crossed=%d rate=%.5f alpha=%.3f verdict=%s\n",
            GUARD_REPLICATIONS, MONITOR_DAYS, guard, guard_rate, ALPHA, if (guard_rate <= ALPHA) "pass" else "fail"))
if (guard_rate > ALPHA) stop("null false-alarm check failed: rate above alpha")

lcg_state <<- LCG_SEED
sticky <- 0L
for (iter in seq_len(STICKY_REPLICATIONS)) {
  control <- 0
  variant <- 0
  crossed_once <- FALSE
  for (day in seq_len(MONITOR_DAYS)) {
    drawn <- draw_sticky_day(0)
    control <- control + drawn[2]
    variant <- variant + drawn[4]
    conversions <- control + variant
    if (conversions > 0 && abs(variant - control) >= mixture_bound(conversions, rhos[reference])) {
      crossed_once <- TRUE
      break
    }
  }
  if (crossed_once) sticky <- sticky + 1L
}
cat(sprintf("metric.check.sticky_assignment block_visitors=%d replications=%d crossed=%d rate=%.5f alpha=%.3f\n",
            STICKY_BLOCK, STICKY_REPLICATIONS, sticky, sticky / STICKY_REPLICATIONS, ALPHA))
lcg_state <<- saved_state  # the example loop resumes the main stream where the invariants left it

for (day in LOOK_POINTS) {
  # %/% binds tighter than * in R: the parentheses keep (250 visitors x day) / 2 per arm.
  n <- (PER_DAY * day) %/% 2L
  cat(sprintf("metric.exact.fixed_date_size horizon_days=%d n_per_arm=%d size=%.6f\n", day, n, exact_size(n)))
}

for (attempt in seq_len(EXAMPLE_LIMIT)) {
  looks <- run_test(PLANNED_DAYS, 0, rhos)
  day <- first_day(abs(looks$z) >= Z_FIXED, PLANNED_DAYS)
  if (day > 0L) {
    cat(sprintf("metric.example.no_effect_declared_by_peeking attempt=%d day=%d z=%.4f gap=%d\n",
                attempt, day, looks$z[day], looks$gap[day]))
    cat(sprintf("metric.example.same_test_at_planned_day day=%d z=%.4f gap=%d significant=%s\n",
                PLANNED_DAYS, looks$z[PLANNED_DAYS], looks$gap[PLANNED_DAYS],
                if (abs(looks$z[PLANNED_DAYS]) >= Z_FIXED) "yes" else "no"))
    cat(sprintf("metric.example.same_test_boundary_crossed_by_planned_day=%s\n",
                if (any(looks$crossed[, reference])) "yes" else "no"))
    break
  }
}

for (scenario in list(list("no_effect", 0), list("true_uplift", UPLIFT))) {
  label <- scenario[[1]]
  uplift <- scenario[[2]]
  fixed <- integer(length(LOOK_POINTS))
  peeks <- integer(length(LOOK_POINTS))
  weekly <- 0L
  after_first_week <- 0L
  crossed_any <- integer(length(rhos))
  crossed_by_day <- matrix(0L, MONITOR_DAYS, length(rhos))
  peek_days <- integer(0)
  boundary_days <- integer(0)
  for (iter in seq_len(REPLICATIONS)) {
    looks <- run_test(MONITOR_DAYS, uplift, rhos)
    significant <- abs(looks$z) >= Z_FIXED
    for (j in seq_along(LOOK_POINTS)) {
      limit <- LOOK_POINTS[j]
      if (significant[limit]) fixed[j] <- fixed[j] + 1L
      day <- first_day(significant, limit)
      if (day > 0L) {
        peeks[j] <- peeks[j] + 1L
        if (limit == PLANNED_DAYS) peek_days <- c(peek_days, day)
      }
    }
    if (any(significant[c(7L, 14L, 21L, PLANNED_DAYS)])) weekly <- weekly + 1L
    if (any(significant[7:PLANNED_DAYS])) after_first_week <- after_first_week + 1L
    for (index in seq_along(rhos)) {
      day <- first_day(looks$crossed[, index], MONITOR_DAYS)
      if (day > 0L) {
        crossed_any[index] <- crossed_any[index] + 1L
        crossed_by_day[day:MONITOR_DAYS, index] <- crossed_by_day[day:MONITOR_DAYS, index] + 1L
        if (index == reference) boundary_days <- c(boundary_days, day)
      }
    }
  }
  for (j in seq_along(LOOK_POINTS)) {
    limit <- LOOK_POINTS[j]
    s <- summary_share(fixed[j], REPLICATIONS)
    cat(sprintf("metric.%s.fixed_date horizon_days=%d declared=%d of %d rate=%.5f mcse=%.5f\n",
                label, limit, s[1], REPLICATIONS, s[2], s[3]))
    s <- summary_share(peeks[j], REPLICATIONS)
    cat(sprintf("metric.%s.daily_peeking horizon_days=%d looks=%d declared=%d of %d rate=%.5f mcse=%.5f\n",
                label, limit, limit, s[1], REPLICATIONS, s[2], s[3]))
  }
  s <- summary_share(weekly, REPLICATIONS)
  cat(sprintf("metric.%s.weekly_reading horizon_days=%d looks=4 declared=%d of %d rate=%.5f mcse=%.5f\n",
              label, PLANNED_DAYS, s[1], REPLICATIONS, s[2], s[3]))
  s <- summary_share(after_first_week, REPLICATIONS)
  cat(sprintf("metric.%s.daily_from_day_7 horizon_days=%d looks=22 declared=%d of %d rate=%.5f mcse=%.5f\n",
              label, PLANNED_DAYS, s[1], REPLICATIONS, s[2], s[3]))
  for (index in seq_along(labels)) {
    s <- summary_share(crossed_any[index], REPLICATIONS)
    cat(sprintf("metric.%s.boundary tuning=%s looks=%d declared=%d of %d rate=%.5f mcse=%.5f\n",
                label, labels[index], MONITOR_DAYS, s[1], REPLICATIONS, s[2], s[3]))
    s <- summary_share(crossed_by_day[PLANNED_DAYS, index], REPLICATIONS)
    cat(sprintf("metric.%s.boundary tuning=%s by_planned_day declared=%d of %d rate=%.5f mcse=%.5f\n",
                label, labels[index], s[1], REPLICATIONS, s[2], s[3]))
  }
  if (label == "true_uplift") {
    goal <- fixed[which(LOOK_POINTS == PLANNED_DAYS)] / REPLICATIONS
    for (index in seq_along(labels)) {
      reached <- 0L
      for (d in seq_len(MONITOR_DAYS)) {
        if (crossed_by_day[d, index] / REPLICATIONS >= goal) { reached <- d; break }
      }
      before_power <- if (reached > 1L) crossed_by_day[reached - 1L, index] / REPLICATIONS else 0
      cat(sprintf("metric.true_uplift.boundary tuning=%s day_matching_fixed_power=%d ratio=%.2f goal=%.5f power_day_before=%.5f\n",
                  labels[index], reached, if (reached > 0L) reached / PLANNED_DAYS else 0, goal, before_power))
    }
    if (length(peek_days)) {
      peek_days <- sort(peek_days)
      cat(sprintf("metric.true_uplift.peeking_median_stop_day=%d within=%d_days\n",
                  peek_days[length(peek_days) %/% 2 + 1], PLANNED_DAYS))
    }
    if (length(boundary_days)) {
      boundary_days <- sort(boundary_days)
      cat(sprintf("metric.true_uplift.boundary_median_stop_day=%d within=%d_days\n",
                  boundary_days[length(boundary_days) %/% 2 + 1], MONITOR_DAYS))
    }
  }
}