SilverSight/julia/PIST/pist_fiedler_chiral.jl

"""
    PIST Fiedler-Aware Chiral Boundary Detection — Julia Port

Extends PIST spectral analysis (SpectralN.lean) with Fiedler vector
sign-pattern analysis for chiral boundary classification.

References:
    - `formal/SilverSight/PIST/SpectralN.lean`
    - `formal/SilverSight/PIST/CartanConnection.lean`
    - `python/pist_fiedler_chiral.py`
"""
module FiedlerChiral

using LinearAlgebra

export build_laplacian_8x8, power_iteration, fiedler_vector,
       classify_chiral_boundary, compute_chiral_boundary_profile

const CHIRAL_LABELS = ["achiral_stable", "left_handed", "right_handed", "chiral_scarred"]

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

function build_laplacian_8x8(cross_coupling::Float64=1e-6)::Matrix{Float64}
    C = zeros(Float64, 8, 8)
    for i in 1:8
        C[i, i] = 39.0 / 256.0
        for j in 1:8
            if i != j
                if div(i - 1, 2) == div(j - 1, 2)
                    C[i, j] = 1.0 / 7.0
                else
                    C[i, j] = cross_coupling
                end
            end
        end
    end

    A = copy(C)
    for i in 1:8; A[i, i] = 0.0; end
    D = diagm(vec(sum(A, dims=2)))
    D - A
end

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

function power_iteration(mat::Matrix{Float64}; max_iter::Int=100, tol::Float64=1e-8)
    n = size(mat, 1)
    v = Float64[Float64(i) for i in 1:n]

    for _ in 1:max_iter
        mv = mat * v
        eig = dot(v, mv) / dot(v, v)
        norm_mv = norm(mv)
        norm_mv < 1e-15 && break
        v_new = mv / norm_mv
        resid = norm(mv - eig * v) / n
        v = v_new
        resid < tol && break
    end

    mv = mat * v
    eig = dot(v, mv) / dot(v, v)
    (eig, v)
end

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

function fiedler_vector(L::Matrix{Float64})
    n = size(L, 1)
    lambda_max, v1 = power_iteration(L)
    # Use full eigendecomposition (n=8 is small enough).
    # For larger n, use iterative methods — for n=8 this is exact.
    eig_vals = eigvals(Symmetric(L))
    eig_vecs = eigvecs(Symmetric(L))

    # Fiedler value = second smallest eigenvalue
    sort_idx = sortperm(eig_vals)
    fiedler_val = eig_vals[sort_idx[2]]
    fiedler_vec = eig_vecs[:, sort_idx[2]]

    (fiedler_val, fiedler_vec)
end

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

function classify_chiral_boundary(fiedler_vec::Vector{Float64})::String
    sign_vec = sign.(fiedler_vec)

    intra_flips = 0
    for k in 0:3
        sign_vec[2k+1] != sign_vec[2k+2] && (intra_flips += 1)
    end

    inter_flips = 0
    for k in 0:2
        sign_vec[2k+2] != sign_vec[2k+3] && (inter_flips += 1)
    end

    bias = [fiedler_vec[2k+1] + fiedler_vec[2k+2] for k in 0:3]
    net_bias = sum(bias)

    if intra_flips == 0 && inter_flips == 0
        return "achiral_stable"
    elseif intra_flips > 0 && net_bias < 0
        return "left_handed"
    elseif intra_flips > 0 && net_bias > 0
        return "right_handed"
    else
        return "chiral_scarred"
    end
end

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

function compute_chiral_boundary_profile(C_matrix::Union{Matrix{Float64}, Nothing}=nothing)
    L = C_matrix === nothing ? build_laplacian_8x8() : build_laplacian_from_matrix(C_matrix)
    f_val, f_vec = fiedler_vector(L)
    chiral_label = classify_chiral_boundary(f_vec)
    lambda_max, _ = power_iteration(L)

    Dict(
        "fiedler_value" => f_val,
        "fiedler_vector" => f_vec,
        "chiral_label" => chiral_label,
        "intra_pair_flips" => sum([sign(f_vec[2k+1]) != sign(f_vec[2k+2]) ? 1 : 0 for k in 0:3]),
        "inter_pair_flips" => sum([sign(f_vec[2k+2]) != sign(f_vec[2k+3]) ? 1 : 0 for k in 0:2]),
        "spectral_gap" => lambda_max - f_val,
        "dominant_eigenvalue" => lambda_max,
    )
end

function build_laplacian_from_matrix(mat::Matrix{Float64})::Matrix{Float64}
    A = abs.(mat)
    for i in 1:size(A, 1); A[i, i] = 0.0; end
    D = diagm(vec(sum(A, dims=2)))
    D - A
end

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

function demo()
    println("="^60)
    println("PIST Fiedler-Aware Chiral Boundary Detection (Julia)")
    println("="^60)

    L = build_laplacian_8x8()
    println("\nLaplacian L:")
    display(round.(L, digits=6))

    lambda_max, v1 = power_iteration(L)
    println("\nλ_max (dominant): $(round(lambda_max, digits=6))")

    f_val, f_vec = fiedler_vector(L)
    println("Fiedler value (λ₂): $(round(f_val, digits=6))")
    println("Spectral gap: $(round(lambda_max - f_val, digits=6))")
    println("Fiedler vector: $(round.(f_vec, digits=6))")
    println("Sign pattern: $(sign.(f_vec))")

    profile = compute_chiral_boundary_profile()
    println("\nChiral classification: $(profile["chiral_label"])")
    println("Intra-pair sign flips: $(profile["intra_pair_flips"])")
    println("Inter-pair sign flips: $(profile["inter_pair_flips"])")

    println("\n--- Perturbation analysis ---")
    L_pert = copy(L)
    L_pert[1, 1] += 0.5
    fv2, fv2_vec = fiedler_vector(L_pert)
    println("Left-bias perturbation: Fiedler=$(round(fv2, digits=6)), chiral=$(classify_chiral_boundary(fv2_vec))")
end

end # module

if abspath(PROGRAM_FILE) == @__FILE__
    using .FiedlerChiral
    FiedlerChiral.demo()
end