|
| 1 | +""" |
| 2 | +scope.py -- Transfer matrix infrastructure for Milestone 6: Scoped Mediation. |
| 3 | +
|
| 4 | +Built by exp_01, imported by exp_02-10. Provides: |
| 5 | +- build_transfer_matrix: construct T mapping parent spectral to child spectral |
| 6 | +- decompose_harmonic_transient: split T = T_harm + T_trans |
| 7 | +- harmonic_fixed_point: iterate T_harm to rank-1 projector |
| 8 | +- scope_attenuation: compute ||T_harm^n|| for n hops |
| 9 | +- pac_budget: compute P, A, xi, Theta from spectral decomposition |
| 10 | +
|
| 11 | +All functions operate on the confluent identity hierarchy built by |
| 12 | +exp_01/exp_02 of the confluent_identity series. |
| 13 | +""" |
| 14 | + |
| 15 | +import numpy as np |
| 16 | +from scipy import sparse |
| 17 | +from scipy.sparse.linalg import eigsh |
| 18 | + |
| 19 | + |
| 20 | +# ============================================================ |
| 21 | +# Constants |
| 22 | +# ============================================================ |
| 23 | +PHI = (1 + np.sqrt(5)) / 2 |
| 24 | +INV_PHI = 1 / PHI |
| 25 | +LN_PHI = np.log(PHI) |
| 26 | +GAMMA_EM = 0.5772156649015329 |
| 27 | +XI_BALANCE = GAMMA_EM + LN_PHI # 1.0584 |
| 28 | + |
| 29 | + |
| 30 | +# ============================================================ |
| 31 | +# Transfer matrix construction |
| 32 | +# ============================================================ |
| 33 | + |
| 34 | +def _get_eigenbasis(L, state_vector, k=10): |
| 35 | + """ |
| 36 | + Compute k eigenvectors of graph Laplacian L, sorted by eigenvalue. |
| 37 | + Returns (eigenvalues, eigenvectors) with zero modes included. |
| 38 | + """ |
| 39 | + n = L.shape[0] |
| 40 | + k_actual = min(k + 1, n - 1) |
| 41 | + |
| 42 | + if k_actual < 2: |
| 43 | + return np.array([0.0]), np.ones((n, 1)) / np.sqrt(n) |
| 44 | + |
| 45 | + if n < 50: |
| 46 | + L_dense = L.toarray() if sparse.issparse(L) else L |
| 47 | + eigenvalues, eigenvectors = np.linalg.eigh(L_dense) |
| 48 | + else: |
| 49 | + try: |
| 50 | + eigenvalues, eigenvectors = eigsh( |
| 51 | + L.astype(float), k=k_actual, which='SM', |
| 52 | + tol=1e-8, maxiter=5000 |
| 53 | + ) |
| 54 | + except Exception: |
| 55 | + L_dense = L.toarray() if sparse.issparse(L) else L |
| 56 | + eigenvalues, eigenvectors = np.linalg.eigh(L_dense) |
| 57 | + eigenvalues = eigenvalues[:k_actual] |
| 58 | + eigenvectors = eigenvectors[:, :k_actual] |
| 59 | + |
| 60 | + idx = np.argsort(eigenvalues) |
| 61 | + return eigenvalues[idx], eigenvectors[:, idx] |
| 62 | + |
| 63 | + |
| 64 | +def build_transfer_matrix(parent_eigvecs, child_indices_in_parent, k=10): |
| 65 | + """ |
| 66 | + Construct the transfer matrix T that maps a parent's spectral identity |
| 67 | + to a child's contribution in that basis. |
| 68 | +
|
| 69 | + T[i,j] = sum_{cell in child} v_i[cell] * v_j[cell] |
| 70 | +
|
| 71 | + where v_i are the parent's eigenvectors. This is the child's "spectral |
| 72 | + footprint" in the parent's eigenbasis -- a k x k matrix whose (i,j) entry |
| 73 | + measures how much the child's cells correlate mode i with mode j. |
| 74 | +
|
| 75 | + Parameters: |
| 76 | + parent_eigvecs: (n_parent, k) eigenvectors of parent's Laplacian |
| 77 | + child_indices_in_parent: local indices of child cells within parent |
| 78 | + k: number of modes to use |
| 79 | +
|
| 80 | + Returns: |
| 81 | + T: (k, k) transfer matrix |
| 82 | + """ |
| 83 | + k_actual = min(k, parent_eigvecs.shape[1]) |
| 84 | + child_vecs = parent_eigvecs[child_indices_in_parent, :k_actual] |
| 85 | + T = child_vecs.T @ child_vecs |
| 86 | + # Normalize by child size so T measures density not total |
| 87 | + T /= max(len(child_indices_in_parent), 1) |
| 88 | + return T |
| 89 | + |
| 90 | + |
| 91 | +def decompose_harmonic_transient(T, n_harmonic=1): |
| 92 | + """ |
| 93 | + Decompose transfer matrix T = T_harm + T_trans. |
| 94 | +
|
| 95 | + T_harm captures the harmonic (zero-mode) component -- the part that |
| 96 | + survives arbitrarily many scope boundaries. T_trans captures the |
| 97 | + transient modes that decay with each hop. |
| 98 | +
|
| 99 | + The harmonic subspace corresponds to the first n_harmonic eigenmodes |
| 100 | + (typically just the zero mode, n_harmonic=1). |
| 101 | +
|
| 102 | + Parameters: |
| 103 | + T: (k, k) transfer matrix |
| 104 | + n_harmonic: number of modes in the harmonic subspace |
| 105 | +
|
| 106 | + Returns: |
| 107 | + T_harm: (k, k) harmonic component |
| 108 | + T_trans: (k, k) transient component |
| 109 | + eigenvalues: eigenvalues of T (sorted descending by magnitude) |
| 110 | + """ |
| 111 | + eigenvalues, eigenvectors = np.linalg.eigh(T) |
| 112 | + # Sort by magnitude (descending) |
| 113 | + idx = np.argsort(np.abs(eigenvalues))[::-1] |
| 114 | + eigenvalues = eigenvalues[idx] |
| 115 | + eigenvectors = eigenvectors[:, idx] |
| 116 | + |
| 117 | + # Harmonic projection: first n_harmonic modes |
| 118 | + V_harm = eigenvectors[:, :n_harmonic] |
| 119 | + T_harm = V_harm @ np.diag(eigenvalues[:n_harmonic]) @ V_harm.T |
| 120 | + |
| 121 | + T_trans = T - T_harm |
| 122 | + return T_harm, T_trans, eigenvalues |
| 123 | + |
| 124 | + |
| 125 | +def harmonic_fixed_point(T_harm, n_iter=20, tol=1e-10): |
| 126 | + """ |
| 127 | + Iterate T_harm^n and test convergence to rank-1 projector. |
| 128 | +
|
| 129 | + Returns: |
| 130 | + converged: bool -- did T_harm^n stabilize? |
| 131 | + rank1_error: float -- ||T^n - T^(n-1)|| at final iteration |
| 132 | + powers: list of (n, T^n, frobenius_norm) at each step |
| 133 | + """ |
| 134 | + T_n = T_harm.copy() |
| 135 | + powers = [(1, T_n.copy(), np.linalg.norm(T_n, 'fro'))] |
| 136 | + |
| 137 | + for n in range(2, n_iter + 1): |
| 138 | + T_prev = T_n.copy() |
| 139 | + T_n = T_n @ T_harm |
| 140 | + norm = np.linalg.norm(T_n, 'fro') |
| 141 | + diff = np.linalg.norm(T_n - T_prev, 'fro') |
| 142 | + powers.append((n, T_n.copy(), norm)) |
| 143 | + |
| 144 | + if diff < tol: |
| 145 | + return True, diff, powers |
| 146 | + |
| 147 | + final_diff = np.linalg.norm(powers[-1][1] - powers[-2][1], 'fro') |
| 148 | + return final_diff < tol, final_diff, powers |
| 149 | + |
| 150 | + |
| 151 | +def scope_attenuation(T_harm, n_hops): |
| 152 | + """ |
| 153 | + Compute the Frobenius norm of T_harm^n for n = 1..n_hops. |
| 154 | +
|
| 155 | + Returns: |
| 156 | + norms: list of float, ||T_harm^n||_F for n=1..n_hops |
| 157 | + ratios: list of float, ||T^(n+1)||/||T^n|| for n=1..n_hops-1 |
| 158 | + """ |
| 159 | + T_n = T_harm.copy() |
| 160 | + norms = [np.linalg.norm(T_n, 'fro')] |
| 161 | + |
| 162 | + for _ in range(n_hops - 1): |
| 163 | + T_n = T_n @ T_harm |
| 164 | + norms.append(np.linalg.norm(T_n, 'fro')) |
| 165 | + |
| 166 | + ratios = [] |
| 167 | + for i in range(len(norms) - 1): |
| 168 | + if norms[i] > 1e-15: |
| 169 | + ratios.append(norms[i + 1] / norms[i]) |
| 170 | + else: |
| 171 | + ratios.append(np.nan) |
| 172 | + |
| 173 | + return norms, ratios |
| 174 | + |
| 175 | + |
| 176 | +def pac_budget(state_vector, L, eigenvectors, eigenvalues): |
| 177 | + """ |
| 178 | + Compute the PAC information budget at a scope boundary. |
| 179 | +
|
| 180 | + P (Potential) = total spectral energy entering |
| 181 | + A (Actualized) = harmonic component (zero-mode projection) |
| 182 | + xi (Structure) = energy in the first few non-zero modes (organized) |
| 183 | + Theta (Thermal) = energy in remaining modes (dissipated) |
| 184 | +
|
| 185 | + Parameters: |
| 186 | + state_vector: field values at cells in this region |
| 187 | + L: graph Laplacian |
| 188 | + eigenvectors: eigenvectors of L |
| 189 | + eigenvalues: eigenvalues of L |
| 190 | +
|
| 191 | + Returns: |
| 192 | + dict with P, A, xi, Theta, conservation_error |
| 193 | + """ |
| 194 | + state_centered = state_vector - np.mean(state_vector) |
| 195 | + coefficients = eigenvectors.T @ state_centered |
| 196 | + energies = coefficients ** 2 |
| 197 | + |
| 198 | + # Total energy |
| 199 | + P = float(np.sum(energies)) |
| 200 | + |
| 201 | + # Harmonic = zero-mode energy |
| 202 | + zero_mask = eigenvalues < 1e-10 |
| 203 | + A = float(np.sum(energies[zero_mask])) |
| 204 | + |
| 205 | + # Structure = energy in first few non-zero modes (modes 1-3) |
| 206 | + nonzero_eigs = np.where(~zero_mask)[0] |
| 207 | + structure_modes = nonzero_eigs[:3] if len(nonzero_eigs) >= 3 else nonzero_eigs |
| 208 | + xi = float(np.sum(energies[structure_modes])) |
| 209 | + |
| 210 | + # Thermal = everything else |
| 211 | + all_budget = set(range(len(eigenvalues))) |
| 212 | + used = set(np.where(zero_mask)[0]) | set(structure_modes) |
| 213 | + thermal_modes = list(all_budget - used) |
| 214 | + Theta = float(np.sum(energies[thermal_modes])) |
| 215 | + |
| 216 | + conservation_error = abs(P - (A + xi + Theta)) |
| 217 | + |
| 218 | + return { |
| 219 | + 'P': P, |
| 220 | + 'A': A, |
| 221 | + 'xi': xi, |
| 222 | + 'Theta': Theta, |
| 223 | + 'conservation_error': conservation_error, |
| 224 | + 'A_fraction': A / P if P > 0 else 0, |
| 225 | + 'xi_fraction': xi / P if P > 0 else 0, |
| 226 | + 'Theta_fraction': Theta / P if P > 0 else 0, |
| 227 | + } |
| 228 | + |
| 229 | + |
| 230 | +def matrix_rank_at_tolerance(M, tol=1e-6): |
| 231 | + """Effective rank of matrix M at given tolerance.""" |
| 232 | + s = np.linalg.svd(M, compute_uv=False) |
| 233 | + return int(np.sum(s > tol * s[0])) |
0 commit comments