1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
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')
|
