#!/usr/bin/env Rscript
#
# NUVMAP Projection Engine — R Reference Implementation
# ======================================================
# Cross-validation: implements the exact same algorithm as the Python
# Q16_16 projection engine, using R's native floating point as a
# third independent reference.  Results are compared against both the
# Python float (Research Stack) and Python Q16_16 (SilverSight) outputs.
#
# Run: Rscript tests/test_nuvmap_r_ref.R

SCALE <- 65536

f2q <- function(f) {
  # Float to Q16_16 raw (external boundary only)
  if (is.nan(f) || is.infinite(f)) return(0L)
  max(-2147483648, min(2147483647, round(f * SCALE)))
}

q2f <- function(q) {
  q / SCALE
}

qdiv <- function(a, b) {
  if (b == 0) return(0L)
  round(as.numeric(a) * SCALE / as.numeric(b))
}

qmul <- function(a, b) {
  round(as.numeric(a) * as.numeric(b) / SCALE)
}

qadd <- function(a, b) as.integer(max(-2147483648, min(2147483647, a + b)))

qsub <- function(a, b) as.integer(max(-2147483648, min(2147483647, a - b)))

# ── Constants ───────────────────────────────────────────────────────

Q16_EPSILON    <- 1L
Q16_ONE        <- 65536L
Q16_HALF       <- 32768L
Q16_ZERO       <- 0L
Q16_PCT1       <- 655L
Q16_PCT70      <- 45875L
Q16_PCT30      <- 19661L
Q16_150PCT     <- 98304L

# ── Test data (same as Python test) ─────────────────────────────────

test_data <- list(
  list(equation_id = 1, amvr = 0.8,  avmr = 0.75, cr = 0.05, cs = "achiral_stable"),
  list(equation_id = 2, amvr = 0.3,  avmr = 0.4,  cr = 0.2,  cs = "left_handed_mass_bias"),
  list(equation_id = 3, amvr = 0.1,  avmr = 0.15, cr = 0.6,  cs = "chiral_scarred"),
  list(equation_id = 4, amvr = 0.5,  avmr = 0.5,  cr = 0.1,  cs = "right_handed_vector_bias"),
  list(equation_id = 5, amvr = 0.9,  avmr = 0.85, cr = 0.02, cs = "achiral_stable")
)

convert_to_q16 <- function(data) {
  lapply(data, function(d) list(
    equation_id = d$equation_id,
    amvr_q16 = f2q(d$amvr),
    avmr_q16 = f2q(d$avmr),
    chiral_residual_q16 = f2q(d$cr),
    chiral_state = d$cs
  ))
}

# ── Projection engine (Python Q16_16 logic) ─────────────────────────

project_q16 <- function(data, total_qubit_budget = 0,
                        chi_max_q16 = Q16_HALF, R_max_q16 = Q16_HALF,
                        landauer_threshold_q16 = Q16_ONE %/% 10) {
  n <- length(data)
  if (n == 0) return(list())

  # max eigenmass
  max_eigenmass <- Q16_ZERO
  for (d in data) {
    raw_q <- qdiv(qadd(d$amvr_q16, d$avmr_q16), Q16_ONE * 2)
    if (raw_q > max_eigenmass) max_eigenmass <- raw_q
  }
  if (max_eigenmass == Q16_ZERO) max_eigenmass <- Q16_ONE

  cells <- list()
  for (i in seq_along(data)) {
    d <- data[[i]]
    amvr_q <- d$amvr_q16
    avmr_q <- d$avmr_q16
    cr_q   <- d$chiral_residual_q16
    eq_id  <- d$equation_id
    cs     <- d$chiral_state

    # E_norm = (amvr + avmr) / 2 / max_eigenmass
    raw_eigenmass <- qdiv(qadd(amvr_q, avmr_q), Q16_ONE * 2)
    E_norm <- qdiv(raw_eigenmass, max_eigenmass)

    # R_i = max(0.01, 1.0 - E_norm)
    one_minus <- qsub(Q16_ONE, E_norm)
    R_i <- max(Q16_PCT1, one_minus)

    # chiral_scarred → 1.5x
    if (cs == "chiral_scarred") {
      R_i <- qmul(R_i, Q16_150PCT)
    }

    # S_i: structural integrity
    if (cs == "achiral_stable") {
      S_i <- Q16_ONE
    } else if (cs %in% c("left_handed_mass_bias", "right_handed_vector_bias")) {
      S_i <- Q16_PCT70
    } else {
      S_i <- Q16_PCT30
    }

    # L_i: Landauer factor
    if (E_norm > landauer_threshold_q16) {
      L_i <- Q16_ONE
    } else {
      L_i <- qdiv(E_norm, landauer_threshold_q16)
    }

    # E_i = lam * v_abs * S_i * L_i / (R_i + epsilon)
    lam_q <- Q16_ONE
    v_abs <- E_norm
    E_i <- qdiv(qmul(qmul(qmul(lam_q, v_abs), S_i), L_i), qadd(R_i, Q16_EPSILON))

    chi_i <- cr_q
    is_R_ok  <- R_i <= R_max_q16
    is_chi_ok <- chi_i <= chi_max_q16
    admissible <- is_R_ok && is_chi_ok

    cells[[i]] <- list(
      u_i = i - 1, v_i = i - 1, k_i = i - 1,
      E_i = E_i, R_i = R_i, chi_i = chi_i,
      S_i = S_i, L_i = L_i, q_i = 0L,
      admissible = admissible,
      equation_id = eq_id
    )
  }

  # Qubit allocation
  total_weight <- Q16_ZERO
  for (c in cells) {
    total_weight <- qadd(total_weight, qdiv(c$E_i, qadd(c$R_i, Q16_EPSILON)))
  }
  if (total_weight == Q16_ZERO) total_weight <- Q16_ONE

  if (total_qubit_budget > 0) {
    budget <- total_qubit_budget
  } else {
    budget <- sum(sapply(cells, function(c) if (c$admissible) floor(c$E_i * 100 / SCALE) else 0))
  }
  budget <- max(budget, sum(sapply(cells, function(c) if (c$admissible) 1 else 0)))

  for (i in seq_along(cells)) {
    if (cells[[i]]$admissible) {
      weight <- qdiv(cells[[i]]$E_i, qadd(cells[[i]]$R_i, Q16_EPSILON))
      raw_q <- floor(budget * weight / total_weight)
      cells[[i]]$q_i <- max(1, raw_q)
    } else {
      cells[[i]]$q_i <- 0L
    }
  }

  # Surface metrics
  total_qubits <- sum(sapply(cells, function(c) c$q_i))
  bekenstein <- if (length(cells) > 0) {
    qdiv(Reduce(`+`, lapply(cells, function(c) c$E_i)), length(cells) * Q16_ONE)
  } else Q16_ZERO
  area_util <- if (bekenstein > 0) {
    qdiv(total_qubits * Q16_ONE, qadd(bekenstein, Q16_EPSILON))
  } else Q16_ZERO

  list(
    cells = cells,
    total_qubits = total_qubits,
    bekenstein_bound = bekenstein,
    area_utilization = area_util
  )
}

# ── Float reference (same logic as Research Stack archive) ──────────

project_float <- function(data, total_qubit_budget = 0,
                          chi_max = 0.5, R_max = 0.5,
                          landauer_threshold = 0.1, eps = 1e-12) {
  n <- length(data)
  if (n == 0) return(list())

  max_eigenmass <- max(sapply(data, function(d) (d$amvr + d$avmr) / 2))
  if (max_eigenmass == 0) max_eigenmass <- 1.0

  cells <- list()
  for (i in seq_along(data)) {
    d <- data[[i]]
    amvr <- d$amvr; avmr <- d$avmr; cr <- d$cr; eq_id <- d$equation_id; cs <- d$cs

    raw_eigenmass <- (amvr + avmr) / 2
    E_norm <- raw_eigenmass / max_eigenmass

    R_i <- max(0.01, 1.0 - E_norm)
    if (cs == "chiral_scarred") R_i <- R_i * 1.5

    S_i <- if (cs == "achiral_stable") 1.0
           else if (cs %in% c("left_handed_mass_bias", "right_handed_vector_bias")) 0.7
           else 0.3

    L_i <- if (E_norm > landauer_threshold) 1.0 else E_norm / landauer_threshold

    lam <- 1.0; v_abs <- E_norm
    E_i <- (lam * v_abs * S_i * L_i) / (R_i + eps)

    chi_i <- cr
    admissible <- (R_i <= R_max) && (chi_i <= chi_max)

    cells[[i]] <- list(
      E_i = E_i, R_i = R_i, chi_i = chi_i,
      S_i = S_i, L_i = L_i, q_i = 0L,
      admissible = admissible,
      equation_id = eq_id
    )
  }

  total_weight <- sum(sapply(cells, function(c) c$E_i / (c$R_i + eps)))
  if (total_weight == 0) total_weight <- 1.0

  if (total_qubit_budget > 0) {
    budget <- total_qubit_budget
  } else {
    budget <- sum(sapply(cells, function(c) if (c$admissible) c$E_i * 100 else 0))
  }
  budget <- max(budget, sum(sapply(cells, function(c) if (c$admissible) 1 else 0)))

  for (i in seq_along(cells)) {
    if (cells[[i]]$admissible) {
      weight <- cells[[i]]$E_i / (cells[[i]]$R_i + eps)
      raw_q <- floor(budget * weight / total_weight)
      cells[[i]]$q_i <- max(1, raw_q)
    } else {
      cells[[i]]$q_i <- 0L
    }
  }

  total_qubits <- sum(sapply(cells, function(c) c$q_i))
  bekenstein <- sum(sapply(cells, function(c) c$E_i)) / max(length(cells), 1)
  area_util <- if (bekenstein > 0) total_qubits / (bekenstein + eps) else 0.0

  list(
    cells = cells,
    total_qubits = total_qubits,
    bekenstein_bound = bekenstein,
    area_utilization = area_util
  )
}

# ── Compare ──────────────────────────────────────────────────────────

compare_all <- function() {
  cat("============================================================\n")
  cat("NUVMAP R Reference — Cross-Validation\n")
  cat("============================================================\n\n")

  data_q16 <- convert_to_q16(test_data)

  # Run R float reference
  cat("R float reference (Research Stack logic):\n")
  rf <- project_float(test_data)
  cat(sprintf("  Cells: %d  Qubits: %d  B: %.6f  U: %.6f\n",
              length(rf$cells), rf$total_qubits, rf$bekenstein_bound, rf$area_utilization))

  # Run R Q16_16 implementation
  cat("\nR Q16_16 (SilverSight logic):\n")
  rq <- project_q16(data_q16)
  cat(sprintf("  Cells: %d  Qubits: %d  B: %.6f  U: %.6f\n",
              length(rq$cells), rq$total_qubits, q2f(rq$bekenstein_bound), q2f(rq$area_utilization)))

  cat("\nCell-by-cell comparison:\n")
  cat(sprintf("  %-5s %-20s %-15s %-15s %-10s %-10s\n",
              "cell", "state", "E_i", "R_i", "q_i", "admit"))
  cat("  ----- -------------------- --------------- --------------- ---------- ----------\n")
  for (i in seq_along(test_data)) {
    fc <- rf$cells[[i]]
    qc <- rq$cells[[i]]
    td <- test_data[[i]]
    cat(sprintf("  [%d]   %-20s float=%-8.4f  q16=%-8.4f  %-4d  %-5s\n",
                i-1, td$cs,
                fc$E_i, q2f(qc$E_i),
                qc$q_i,
                if (qc$admissible) "OK" else "REJ"))
    cat(sprintf("        %-20s R=%-8.6f      R=%-8.6f\n",
                "", fc$R_i, q2f(qc$R_i)))
  }

  # Summary check
  cat("\n--- Tolerance check (< 1% relative divergence except ε effect) ---\n")
  all_ok <- TRUE
  for (i in seq_along(test_data)) {
    fc <- rf$cells[[i]]
    qc <- rq$cells[[i]]
    e_diff <- abs(q2f(qc$E_i) - fc$E_i)
    e_rel  <- e_diff / max(abs(fc$E_i), 1e-12)
    if (e_rel > 0.02) {
      cat(sprintf("  ✗ cell[%d] E_i rel diff: %.4f (%.6f vs %.6f)\n",
                  i-1, e_rel, fc$E_i, q2f(qc$E_i)))
      all_ok <- FALSE
    }
    q_diff <- abs(qc$q_i - fc$q_i)
    q_rel  <- q_diff / max(fc$q_i, 1)
    if (q_diff > 1 && q_rel > 0.01) {
      cat(sprintf("  ✗ cell[%d] q_i: %d vs %d (rel=%.4f)\n",
                  i-1, fc$q_i, qc$q_i, q_rel))
      all_ok <- FALSE
    }
  }
  if (all_ok) cat("  ✓ All cells within tolerance\n")

  cat("\n============================================================\n")
  if (all_ok) {
    cat("R REFERENCE: PASS — Q16_16 implementation verified\n")
    quit(save = "no", status = 0)
  } else {
    cat("R REFERENCE: FAIL — numerical divergence detected\n")
    quit(save = "no", status = 1)
  }
}

compare_all()