#!/usr/bin/env Rscript
# MSC-P-033 marketing forecast validation. MIT License. Base R reference.
# The design matrix, the origins, the horizon, the nominal level and both
# thresholds are declared, so this file reproduces the same digits as the
# Python reference. Nothing is random.

args <- commandArgs(trailingOnly = TRUE)
path <- if (length(args) >= 1) args[[1]] else file.path("..", "datasets", "msc-p033-weekly-series.csv")

raw_lines <- readLines(path, warn = FALSE)
raw_lines <- raw_lines[nzchar(raw_lines)]
if (raw_lines[[1]] != "week,revenue_eur") stop("exact schema required")
if (any(lengths(regmatches(raw_lines[-1], gregexpr(",", raw_lines[-1], fixed = TRUE))) != 1L)) stop("exactly two cells required per row")

d <- read.csv(path, stringsAsFactors = FALSE, colClasses = c("integer", "numeric"))
if (!identical(names(d), c("week", "revenue_eur"))) stop("exact schema required")

PERIOD <- 52L
HARMONICS <- 2L
FIRST_ORIGIN <- 104L
ORIGIN_STEP <- 4L
LAST_ORIGIN <- 152L
HORIZON <- 4L
NOMINAL <- 0.80
NORMAL_QUANTILE <- 1.281552  # the 0.90 quantile of the standard normal law
MASE_MAX <- 1.00
COVERAGE_TOLERANCE <- 0.10

if (nrow(d) < LAST_ORIGIN + HORIZON) stop("the series is shorter than the declared validation design")
if (!identical(d$week, seq_len(nrow(d)))) stop("weeks must be numbered 1, 2, ... without gaps")
if (any(!is.finite(d$revenue_eur)) || any(d$revenue_eur <= 0)) stop("the outcome must be finite and strictly positive")
y <- d$revenue_eur

design <- function(weeks) {
  columns <- list(rep(1, length(weeks)), as.numeric(weeks))
  for (k in seq_len(HARMONICS)) {
    angle <- 2 * pi * k * weeks / PERIOD
    columns <- c(columns, list(sin(angle)), list(cos(angle)))
  }
  do.call(cbind, columns)
}

interval_score <- function(lower, upper, actual, alpha) {
  score <- upper - lower
  if (actual < lower) score <- score + 2 / alpha * (lower - actual)
  if (actual > upper) score <- score + 2 / alpha * (actual - upper)
  score
}

alpha <- 1 - NOMINAL
origins <- seq(FIRST_ORIGIN, LAST_ORIGIN, by = ORIGIN_STEP)
scaled_model <- c()
scaled_benchmark <- c()
per_horizon <- vector("list", HORIZON)
for (h in seq_len(HORIZON)) per_horizon[[h]] <- numeric(0)
covered <- 0L
widths <- c()
scores <- c()
absolute <- c()

for (origin in origins) {
  weeks <- seq_len(origin)
  X <- design(weeks)
  normal <- crossprod(X)
  if (rcond(normal) < 1e-12) stop("the declared design matrix is singular on this window")
  coefficients <- as.vector(solve(normal, crossprod(X, y[weeks])))
  residuals <- y[weeks] - as.vector(X %*% coefficients)
  degrees <- origin - ncol(X)
  if (degrees <= 0) stop("too few observations before the origin for the declared design")
  deviation <- sqrt(sum(residuals^2) / degrees)
  seasonal_errors <- abs(y[(PERIOD + 1):origin] - y[1:(origin - PERIOD)])
  scale_value <- mean(seasonal_errors)
  if (!is.finite(scale_value) || scale_value <= 0) stop("the declared seasonal scale is zero")
  half_width <- NORMAL_QUANTILE * deviation
  for (h in seq_len(HORIZON)) {
    week <- origin + h
    actual <- y[week]
    prediction <- as.vector(design(week) %*% coefficients)
    error <- abs(actual - prediction)
    absolute <- c(absolute, error)
    scaled_model <- c(scaled_model, error / scale_value)
    per_horizon[[h]] <- c(per_horizon[[h]], error / scale_value)
    scaled_benchmark <- c(scaled_benchmark, abs(actual - y[week - PERIOD]) / scale_value)
    lower <- prediction - half_width
    upper <- prediction + half_width
    if (lower <= actual && actual <= upper) covered <- covered + 1L
    widths <- c(widths, upper - lower)
    scores <- c(scores, interval_score(lower, upper, actual, alpha))
  }
}

count <- length(scaled_model)
mase <- mean(scaled_model)
benchmark <- mean(scaled_benchmark)
coverage <- covered / count
accuracy_flag <- if (mase < MASE_MAX) "PASS" else "FAIL"
calibration_flag <- if (abs(coverage - NOMINAL) <= COVERAGE_TOLERANCE) "PASS" else "FAIL"
verdict <- if (accuracy_flag == "PASS" && calibration_flag == "PASS") "FORECAST_READABLE_FOR_PLANNING" else "DIAGNOSTIC_BLOCKS_FORECAST_READING"

cat(sprintf("origins=%d\n", length(origins)))
cat(sprintf("first_origin=%d\n", origins[[1]]))
cat(sprintf("last_origin=%d\n", origins[[length(origins)]]))
cat(sprintf("horizon=%d\n", HORIZON))
cat(sprintf("forecasts=%d\n", count))
cat(sprintf("model_mase=%.6f\n", mase))
cat(sprintf("benchmark_mase=%.6f\n", benchmark))
cat(sprintf("model_mae=%.6f\n", mean(absolute)))
cat(sprintf("accuracy_flag=%s\n", accuracy_flag))
for (h in seq_len(HORIZON)) cat(sprintf("mase_h%d=%.6f\n", h, mean(per_horizon[[h]])))
cat(sprintf("nominal_coverage=%.2f\n", NOMINAL))
cat(sprintf("empirical_coverage=%.6f\n", coverage))
cat(sprintf("mean_interval_width=%.6f\n", mean(widths)))
cat(sprintf("mean_interval_score=%.6f\n", mean(scores)))
cat(sprintf("calibration_flag=%s\n", calibration_flag))
cat(sprintf("verdict=%s\n", verdict))