# PIST Fiedler-Aware Chiral Boundary Detection — R Port
#
# Mirrors `python/pist_fiedler_chiral.py` and `julia/PIST/pist_fiedler_chiral.jl`.
#
# Extends PIST spectral analysis (SpectralN.lean) with Fiedler vector
# sign-pattern analysis for chiral boundary classification.

# ── Build Laplacian ──────────────────────────────────────────────────

build_laplacian_8x8 <- function(cross_coupling = 1e-6) {
  C <- matrix(0, nrow = 8, ncol = 8)
  for (i in 1:8) {
    C[i, i] <- 39 / 256
    for (j in 1:8) {
      if (i != j) {
        if (floor((i - 1) / 2) == floor((j - 1) / 2)) {
          C[i, j] <- 1 / 7
        } else {
          C[i, j] <- cross_coupling
        }
      }
    }
  }

  A <- C
  diag(A) <- 0
  D <- diag(rowSums(A))
  D - A
}

# ── Power Iteration ──────────────────────────────────────────────────

power_iteration <- function(mat, max_iter = 100, tol = 1e-8) {
  n <- nrow(mat)
  v <- as.numeric(1:n)

  for (iter in 1:max_iter) {
    mv <- mat %*% v
    eig <- as.numeric((t(v) %*% mv) / (t(v) %*% v))
    norm_mv <- sqrt(sum(mv^2))
    if (norm_mv < 1e-15) break
    v_new <- as.numeric(mv / norm_mv)
    resid <- sqrt(sum((mv - eig * v)^2)) / n
    v <- v_new
    if (resid < tol) break
  }

  mv <- mat %*% v
  eig <- as.numeric((t(v) %*% mv) / (t(v) %*% v))
  list(eigenvalue = eig, eigenvector = as.numeric(v))
}

# ── Fiedler Vector ───────────────────────────────────────────────────

fiedler_vector <- function(L) {
  n <- nrow(L)
  res <- power_iteration(L)
  lambda_max <- res$eigenvalue
  mu <- lambda_max
  shift_mat <- mu * diag(n) - L

  v <- rep(1, n) / sqrt(n)

  for (iter in 1:50) {
    w <- tryCatch(solve(shift_mat, v), error = function(e) NULL)
    if (is.null(w)) break
    norm_w <- sqrt(sum(w^2))
    if (norm_w < 1e-10) break
    v_new <- as.numeric(w / norm_w)
    eig <- as.numeric(t(v_new) %*% (L %*% v_new))
    if (sqrt(sum((v_new - v)^2)) < 1e-8) {
      v <- v_new
      break
    }
    v <- v_new
  }

  fiedler_val <- as.numeric(t(v) %*% (L %*% v))
  list(fiedler_value = fiedler_val, fiedler_vector = v)
}

# ── Chiral Classification ────────────────────────────────────────────

classify_chiral_boundary <- function(fiedler_vec) {
  s <- sign(fiedler_vec)

  intra_flips <- 0
  for (k in 0:3) {
    if (s[2*k + 1] != s[2*k + 2]) intra_flips <- intra_flips + 1
  }

  inter_flips <- 0
  for (k in 0:2) {
    if (s[2*k + 2] != s[2*k + 3]) inter_flips <- inter_flips + 1
  }

  bias <- sapply(0:3, function(k) fiedler_vec[2*k + 1] + fiedler_vec[2*k + 2])
  net_bias <- sum(bias)

  if (intra_flips == 0 && inter_flips == 0) {
    return("achiral_stable")
  } else if (intra_flips > 0 && net_bias < 0) {
    return("left_handed")
  } else if (intra_flips > 0 && net_bias > 0) {
    return("right_handed")
  } else {
    return("chiral_scarred")
  }
}

# ── Full Profile ─────────────────────────────────────────────────────

compute_chiral_boundary_profile <- function(C_matrix = NULL) {
  L <- if (is.null(C_matrix)) build_laplacian_8x8() else build_laplacian_from_matrix(C_matrix)
  fres <- fiedler_vector(L)
  chiral_label <- classify_chiral_boundary(fres$fiedler_vector)
  pres <- power_iteration(L)
  s <- sign(fres$fiedler_vector)

  list(
    fiedler_value = fres$fiedler_value,
    fiedler_vector = fres$fiedler_vector,
    chiral_label = chiral_label,
    intra_pair_flips = sum(sapply(0:3, function(k) if (s[2*k+1] != s[2*k+2]) 1 else 0)),
    inter_pair_flips = sum(sapply(0:2, function(k) if (s[2*k+2] != s[2*k+3]) 1 else 0)),
    spectral_gap = pres$eigenvalue - fres$fiedler_value,
    dominant_eigenvalue = pres$eigenvalue
  )
}

build_laplacian_from_matrix <- function(mat) {
  A <- abs(mat)
  diag(A) <- 0
  D <- diag(rowSums(A))
  D - A
}

# ── Demo ──────────────────────────────────────────────────────────────

demo <- function() {
  cat(paste(rep("=", 60), collapse = ""), "\n")
  cat("PIST Fiedler-Aware Chiral Boundary Detection (R)\n")
  cat(paste(rep("=", 60), collapse = ""), "\n")

  L <- build_laplacian_8x8()
  cat("\nLaplacian L:\n")
  print(round(L, 6))

  pres <- power_iteration(L)
  cat(sprintf("\nλ_max (dominant): %.6f\n", pres$eigenvalue))

  fres <- fiedler_vector(L)
  cat(sprintf("Fiedler value (λ₂): %.6f\n", fres$fiedler_value))
  cat(sprintf("Spectral gap: %.6f\n", pres$eigenvalue - fres$fiedler_value))
  cat("Fiedler vector:", sprintf("%.6f", fres$fiedler_vector), "\n")
  cat("Sign pattern:", sprintf("%+d", sign(fres$fiedler_vector)), "\n")

  profile <- compute_chiral_boundary_profile()
  cat(sprintf("\nChiral classification: %s\n", profile$chiral_label))
  cat(sprintf("Intra-pair sign flips: %d\n", profile$intra_pair_flips))
  cat(sprintf("Inter-pair sign flips: %d\n", profile$inter_pair_flips))

  cat("\n--- Perturbation analysis ---\n")
  L_pert <- L
  L_pert[1, 1] <- L_pert[1, 1] + 0.5
  fv2 <- fiedler_vector(L_pert)
  cat(sprintf("Left-bias perturbation: Fiedler=%.6f, chiral=%s\n",
      fv2$fiedler_value, classify_chiral_boundary(fv2$fiedler_vector)))
}

if (interactive() && Sys.getenv("R_TEST") == "") {
  demo()
}