SilverSight/tests/test_nuvmap_r_ref.R

327 lines
10 KiB
R

#!/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()