diff options
Diffstat (limited to 'hw_sim/toy_loop.py')
| -rw-r--r-- | hw_sim/toy_loop.py | 160 |
1 files changed, 160 insertions, 0 deletions
diff --git a/hw_sim/toy_loop.py b/hw_sim/toy_loop.py new file mode 100644 index 0000000..6e4e7c9 --- /dev/null +++ b/hw_sim/toy_loop.py @@ -0,0 +1,160 @@ +"""Stage C: SPICE-in-the-loop training of a toy 2-layer net (8 -> 8 tanh -> 4 softmax). +Every MVM — forward AND transpose (the error transport J^T read) — is solved by ngspice on a +differential resistive column bank with: per-bit ladder mismatch (fixed per device), per-column +input offsets, 8-bit code quantization of a digital fp32 master (the T64 word-streaming story), +and the v1.1-measured read noise (139 uV) added at each read. +Three arms: SPICE loop | ideal numpy | behavioral numpy (same non-idealities, algebraic). +Curve overlap SPICE ~= behavioral certifies the fidelity ladder used at 72M+. +""" +import os, time +os.environ.setdefault('NGSPICE_LIBRARY_PATH', '/home/yurenh2/miniconda3/lib/libngspice.so') +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from PySpice.Spice.Netlist import Circuit + +rng = np.random.default_rng(3) +D_IN, D_H, D_OUT = 8, 8, 4 +R0 = 11e3 +RF = 1e3 +VOS_SIG = 500e-6 +READ_NOISE = 139e-6 +SIG_BIT = 0.004 +STEPS = 300 +BATCH = 8 +LR = 0.15 +BETA_SCALE = 1.0 # top nudge folded into d2 analytically (single-layer-linear blocks -> K=1 exact) + +# ---- the "device": fixed mismatch/offset samples ---- +def make_device(rows, cols): + return { + 'eps_p': rng.normal(0, SIG_BIT, (rows, cols, 8)), + 'eps_n': rng.normal(0, SIG_BIT, (rows, cols, 8)), + 'vos_col': rng.normal(0, VOS_SIG, cols), + 'vos_row': rng.normal(0, VOS_SIG, rows), + } + +DEV1 = make_device(D_IN, D_H) +DEV2 = make_device(D_H, D_OUT) + +def code_eff(codes, eps): + bits = ((codes[..., None].astype(int) >> np.arange(8)) & 1) + w = (2 ** np.arange(8)) * (1 + eps) + return (bits * w).sum(-1) / 256.0 + +def quant_codes(W): + Wp = np.clip(W, 0, None); Wn = np.clip(-W, 0, None) + scale = max(np.abs(W).max(), 1e-6) + return (np.round(Wp / scale * 255).astype(int), np.round(Wn / scale * 255).astype(int), scale) + +def spice_mvm(W, x, dev, transpose=False): + """y = W^T x (forward, x over rows) or y = W d (transpose, d over cols) via ngspice DC.""" + cp, cn, scale = quant_codes(W) + gp = code_eff(cp, dev['eps_p']) / R0 # rows x cols conductances + gn = code_eff(cn, dev['eps_n']) / R0 + rows, cols = W.shape + c = Circuit('mvm') + if not transpose: + for i in range(rows): + c.V(f'x{i}', f'r{i}', c.gnd, float(x[i])) + outs = [] + for j in range(cols): + for i in range(rows): + if gp[i, j] > 1e-9: c.R(f'p{i}_{j}', f'r{i}', f'sp{j}', 1.0 / gp[i, j]) + if gn[i, j] > 1e-9: c.R(f'n{i}_{j}', f'r{i}', f'sn{j}', 1.0 / gn[i, j]) + for tag in ('sp', 'sn'): + c.B(f'amp{tag}{j}', f'o{tag}{j}', c.gnd, v=f'-1e6*(v({tag}{j})-{dev["vos_col"][j]})') + c.R(f'f{tag}{j}', f'o{tag}{j}', f'{tag}{j}', RF) + outs.append(j) + sim = c.simulator(temperature=27, nominal_temperature=27) + an = sim.operating_point() + y = np.array([float(an[f'osp{j}']) - float(an[f'osn{j}']) for j in outs]) + else: + d = x + for j in range(cols): + c.V(f'd{j}', f'c{j}', c.gnd, float(d[j])) + for i in range(rows): + for j in range(cols): + if gp[i, j] > 1e-9: c.R(f'p{i}_{j}', f'c{j}', f'sp{i}', 1.0 / gp[i, j]) + if gn[i, j] > 1e-9: c.R(f'n{i}_{j}', f'c{j}', f'sn{i}', 1.0 / gn[i, j]) + for i in range(rows): + for tag in ('sp', 'sn'): + c.B(f'amp{tag}{i}', f'o{tag}{i}', c.gnd, v=f'-1e6*(v({tag}{i})-{dev["vos_row"][i]})') + c.R(f'f{tag}{i}', f'o{tag}{i}', f'{tag}{i}', RF) + sim = c.simulator(temperature=27, nominal_temperature=27) + an = sim.operating_point() + y = np.array([float(an[f'osp{i}']) - float(an[f'osn{i}']) for i in range(rows)]) + y = -y / RF * R0 * scale # undo -RF*G scaling back to W-units + return y + rng.normal(0, READ_NOISE / RF * R0 * scale, y.shape) + +def behav_mvm(W, x, dev, transpose=False): + cp, cn, scale = quant_codes(W) + gp = code_eff(cp, dev['eps_p']); gn = code_eff(cn, dev['eps_n']) + Weff = (gp - gn) * scale + y = (Weff.T @ x) if not transpose else (Weff @ x) + return y + rng.normal(0, READ_NOISE / RF * R0 * scale, y.shape) + +def ideal_mvm(W, x, dev, transpose=False): + return (W.T @ x) if not transpose else (W @ x) + +# ---- data: 4 Gaussian blobs in 8-D ---- +NTR = 256 +centers = rng.normal(0, 1.0, (D_OUT, D_IN)) +Xtr = np.concatenate([c + 0.35 * rng.normal(0, 1, (NTR // D_OUT, D_IN)) for c in centers]) +Ytr = np.concatenate([np.full(NTR // D_OUT, k) for k in range(D_OUT)]) +perm = rng.permutation(NTR); Xtr, Ytr = Xtr[perm], Ytr[perm] + +def softmax(z): + e = np.exp(z - z.max()); return e / e.sum() + +def train(mvm, tag): + r = np.random.default_rng(42) + W1 = r.normal(0, 0.3, (D_IN, D_H)) + W2 = r.normal(0, 0.3, (D_H, D_OUT)) + losses = [] + t0 = time.time() + for step in range(STEPS): + idx = r.integers(0, NTR, BATCH) + gW1 = np.zeros_like(W1); gW2 = np.zeros_like(W2); L = 0.0 + for s in idx: + x, yl = Xtr[s], Ytr[s] + a1 = mvm(W1, x, DEV1) # analog MVM 1 + h = np.tanh(a1) # digital activation (residency theorem) + a2 = mvm(W2, h, DEV2) # analog MVM 2 + p = softmax(a2) + L -= np.log(max(p[yl], 1e-9)) + d2 = p.copy(); d2[yl] -= 1.0 # top nudge (digital loss head) + d1 = mvm(W2, d2, DEV2, transpose=True) * (1 - h ** 2) # ANALOG transpose read + gW2 += np.outer(h, d2); gW1 += np.outer(x, d1) + W1 -= LR * gW1 / BATCH; W2 -= LR * gW2 / BATCH + losses.append(L / BATCH) + if step % 50 == 0: + print(f'[{tag}] step {step} loss {losses[-1]:.4f} ({time.time()-t0:.0f}s)', flush=True) + # final train accuracy + acc = 0 + for s in range(NTR): + h = np.tanh(mvm(W1, Xtr[s], DEV1)); p = softmax(mvm(W2, h, DEV2)) + acc += int(p.argmax() == Ytr[s]) + print(f'[{tag}] final loss {np.mean(losses[-20:]):.4f} | train acc {acc/NTR*100:.1f}% ' + f'| wall {time.time()-t0:.0f}s', flush=True) + return np.array(losses), acc / NTR + +l_ideal, a_ideal = train(ideal_mvm, 'ideal') +l_behav, a_behav = train(behav_mvm, 'behavioral') +l_spice, a_spice = train(spice_mvm, 'SPICE') + +fig, ax = plt.subplots(figsize=(8.6, 4.8)) +k = np.ones(10) / 10 +for l, name, col in ((l_ideal, f'ideal numpy (acc {a_ideal*100:.0f}%)', '#888888'), + (l_behav, f'behavioral non-idealities (acc {a_behav*100:.0f}%)', '#2c6fbb'), + (l_spice, f'SPICE-in-the-loop (acc {a_spice*100:.0f}%)', '#d95f02')): + ax.plot(np.convolve(l, k, 'valid'), label=name, lw=1.6, + color=col, alpha=0.9) +ax.set_xlabel('training step'); ax.set_ylabel('train CE (smoothed)') +ax.set_title('Stage C: every MVM (forward + transpose) through ngspice — 8-bit codes, mismatch,\n' + 'offsets, measured read noise. Ladder certified if SPICE ≈ behavioral.') +ax.legend(fontsize=9) +fig.tight_layout() +fig.savefig('/home/yurenh2/ept/assets/figs/fig_spice_loop.png', dpi=150) +print('DONE_SPICE_LOOP') |