1import json
  2import numpy as np
  3from scipy.linalg import expm, eigvalsh
  4from scipy.optimize import linear_sum_assignment
  5
  6SEED = 1502
  7rng = np.random.default_rng(SEED)
  8
  9
 10def reachable(MA, MB):
 11    n = MA.shape[0]
 12    seen = np.any(MB != 0, axis=1)
 13    changed = True
 14    while changed:
 15        changed = False
 16        newly = (MA.astype(bool) @ seen.astype(int)) > 0
 17        add = newly & ~seen
 18        if np.any(add): seen |= add; changed = True
 19    return seen
 20
 21
 22def boolean_core(MA, MB):
 23    """Boolean sparsity of [B, AB, ..., A^(n-1)B], retaining each block."""
 24    n, m = MB.shape
 25    P = MB.astype(bool)
 26    blocks = [P.copy()]
 27    for _ in range(1, n):
 28        P = (MA.astype(int) @ P.astype(int)) > 0
 29        blocks.append(P.copy())
 30    return np.concatenate(blocks, axis=1)
 31
 32
 33def max_matching_rows(pattern):
 34    r, c = pattern.shape
 35    if not np.any(pattern): return 0, []
 36    # Add a large penalty for forbidden entries; assignment then maximizes allowed pairs.
 37    cost = np.where(pattern, 0.0, 1e6)
 38    ri, ci = linear_sum_assignment(cost)
 39    pairs = [(int(i), int(j)) for i, j in zip(ri, ci) if pattern[i, j]]
 40    return len(pairs), pairs
 41
 42
 43def ctrb(A, B):
 44    n = A.shape[0]; P = B.copy(); blocks = []
 45    for _ in range(n): blocks.append(P); P = A @ P
 46    return np.concatenate(blocks, axis=1)
 47
 48
 49def gramian(A, B, T=8., steps=240):
 50    dt = T / steps; W = np.zeros((A.shape[0], A.shape[0]))
 51    for q in range(steps):
 52        E = expm(A * ((q + .5) * dt)); X = E @ B
 53        W += dt * (X @ X.T)
 54    return (W + W.T) / 2
 55
 56
 57def normalized_design(kind, n=10, m=2):
 58    MA = np.zeros((n, n), int); MB = np.zeros((n, m), int)
 59    if kind == 'matched':
 60        MB[0, 0] = 1; MB[5, 1] = 1
 61        for i in range(1, 5): MA[i, i-1] = 1
 62        for i in range(6, 10): MA[i, i-1] = 1
 63    elif kind == 'unmatched':
 64        # Every state is reachable, but all descendants have the same predecessor/input history.
 65        MB[0, 0] = 1
 66        for i in range(1, n): MA[i, 0] = 1
 67    elif kind == 'inaccessible':
 68        MB[0, 0] = 1
 69        for i in range(1, n-1): MA[i, i-1] = 1
 70        MA[n-1, n-1] = 1
 71    elif kind == 'random':
 72        MA = (rng.random((n, n)) < .22).astype(int); np.fill_diagonal(MA, 0)
 73        MB = (rng.random((n, m)) < .5).astype(int)
 74        if not np.any(MB): MB[0, 0] = 1
 75    return MA, MB
 76
 77
 78def weighted(MA, MB, seed, stable=True):
 79    r = np.random.default_rng(seed)
 80    A = MA * r.normal(0, .35, MA.shape); B = MB * r.normal(0, 1, MB.shape)
 81    rad = max(abs(np.linalg.eigvals(A))) if np.any(A) else 0
 82    if stable and rad > 0: A *= .65 / rad
 83    return A, B
 84
 85
 86def rank_tol(M):
 87    s = np.linalg.svd(M, compute_uv=False)
 88    return int(np.sum(s > 1e-9 * max(s[0], 1e-30)))
 89
 90
 91def structural_report():
 92    out = []
 93    for kind in ['matched', 'unmatched', 'inaccessible', 'random']:
 94        MA, MB = normalized_design(kind); core = boolean_core(MA, MB); mm, _ = max_matching_rows(core)
 95        out.append({'design': kind, 'reachable': int(reachable(MA, MB).sum()), 'n': len(MA),
 96                    'core_nonzero': int(core.sum()), 'matching': mm, 'row_saturating': mm == len(MA)})
 97    return out
 98
 99
100def gramian_sweep():
101    rows = []
102    for kind in ['matched', 'unmatched', 'inaccessible', 'random']:
103        MA, MB = normalized_design(kind); ranks = []; mins = []; ctranks = []
104        for seed in range(8):
105            A, B = weighted(MA, MB, seed); ctranks.append(rank_tol(ctrb(A, B)))
106            W = gramian(A, B, T=8., steps=180); ev = eigvalsh(W)
107            ranks.append(int(np.sum(ev > 1e-8 * max(ev[-1], 1e-30)))); mins.append(float(ev[0]))
108        rows.append({'design': kind, 'mean_ctrb_rank': float(np.mean(ctranks)),
109                     'mean_gramian_rank': float(np.mean(ranks)), 'min_gramian_rank': min(ranks),
110                     'mean_lambda_min': float(np.mean(mins)), 'median_lambda_min': float(np.median(mins))})
111    return rows
112
113
114def horizon_sweep():
115    MA, MB = normalized_design('matched'); A, B = weighted(MA, MB, 7); result = []
116    for T in [.25, .5, 1, 2, 4, 8]:
117        W = gramian(A, B, T=T, steps=180); ev = eigvalsh(W)
118        result.append({'T': T, 'lambda_min': float(ev[0]), 'lambda_max': float(ev[-1]),
119                       'normalized_lambda_min': float(ev[0] / ev[-1]),
120                       'rank': int(np.sum(ev > 1e-8 * ev[-1]))})
121    return result
122
123
124
125def horizon_controllability_sweep():
126    # Discrete-time analogue: rank grows when additional input history reaches
127    # successive chain states; the unmatched star saturates at its input rank.
128    rows=[]
129    for kind in ["matched", "unmatched", "inaccessible"]:
130        MA,MB=normalized_design(kind); A,B=weighted(MA,MB,7)
131        for H in [1,2,3,4,5,6,8,10]:
132            P=B.copy(); blocks=[]
133            for _ in range(H): blocks.append(P); P=A@P
134            C=np.concatenate(blocks,axis=1)
135            sv=np.linalg.svd(C,compute_uv=False)
136            rows.append({"design":kind,"horizon":H,"rank":rank_tol(C),
137                         "min_nonzero_singular":float(sv[rank_tol(C)-1]) if rank_tol(C)>0 else 0.0})
138    return rows
139
140
141def baseline_comparison():
142    # Standard sparse random initialization versus matching mask, same n,m and
143    # one fixed linear recurrent setup; report the controllability metric.
144    MA,MB=normalized_design("matched"); RA,RB=normalized_design("random")
145    vals=[]
146    for seed in range(8):
147        for name,ma,mb in [("matching",MA,MB),("random_sparse",RA,RB)]:
148            A,B=weighted(ma,mb,seed); W=gramian(A,B,T=8.,steps=180); ev=eigvalsh(W)
149            vals.append({"design":name,"rank":int(np.sum(ev>1e-8*max(ev[-1],1e-30))),
150                         "lambda_min":float(max(ev[0],0.0))})
151    return vals
152
153def main():
154    result = {'seed': SEED, 'structural': structural_report(), 'gramian_sweep': gramian_sweep(),
155              'matched_horizon_sweep': horizon_sweep(),
156              'horizon_controllability_sweep': horizon_controllability_sweep(),
157              'baseline_comparison': baseline_comparison(),
158              'predictions': [
159                'accessibility plus row-saturating matching predicts generic full controllability rank',
160                'the matched finite-horizon Gramian has positive minimum eigenvalue for generic weights',
161                'increasing horizon increases the minimum Gramian eigenvalue for the matched chain, while inaccessible remains singular']}
162    with open('results.json', 'w') as f: json.dump(result, f, indent=2)
163    print(json.dumps(result, indent=2))
164
165if __name__ == '__main__': main()