Shin et al. (2025) from dandiset 000248: checks of the decoder and other decoder types¶

Dandiset 000248 version 0.240502.2344. This notebook repeats decoder reading (a) of decoding_shin_2025.ipynb with other cross-validation seeds (for fig3c-npx-v1-inference and fig3d-npx-areas-inference; reading (b) with the same seeds is in decoding_seed1_shin_2025.ipynb and decoding_seed2_shin_2025.ipynb), checks the bound of the C search and the coding of decoder reading (b) in V1, PM and AM, and tests supp3-decoder-types.

Claim Tested here
fig3c-npx-v1-inference, fig3d-npx-areas-inference seeds of reading (a) only; outcomes in decoding_shin_2025.ipynb
supp3-decoder-types quadratic and RBF SVMs; the network is timed only
In [1]:
import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import time
import warnings

import numpy as np
import pandas as pd
from joblib import Parallel, delayed
from scipy import stats

import illusion as il

warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
N_JOBS = 4

started = time.time()
sessions = il.load_all(ccg=False)
print(f"loaded {len(sessions)} sessions of dandiset {il.DANDISET_ID} version {il.DANDISET_VERSION} "
      f"in {time.time() - started:.0f} s")
loaded 12 sessions of dandiset 000248 version 0.240502.2344 in 133 s

Data check¶

RS units per area and session, and ICwcfg1 trials of the decoded images.

In [2]:
rows = []
for s in sessions:
    u, b = s["units"], s["blocks"]["ICwcfg1"]
    n = {c: int(np.sum(b["codes"] == c)) for c in [il.IC1, il.LC1, il.LC2, il.IC2, il.TRE1, il.TRE2, il.XRE1, il.XRE2]}
    rows.append({"subject": s["subject"], **{f"{a}_RS": int(((u.area == a) & u.rs).sum()) for a in il.AREAS},
                 "trials IC1/LC1/LC2/IC2": "/".join(str(n[c]) for c in [il.IC1, il.LC1, il.LC2, il.IC2]),
                 "TRE1/TRE2": f"{n[il.TRE1]}/{n[il.TRE2]}", "XRE1/XRE2": f"{n[il.XRE1]}/{n[il.XRE2]}"})
units_tab = pd.DataFrame(rows)
print(units_tab.to_string(index=False))
allR = np.concatenate([s["blocks"]["ICwcfg1"]["R"].ravel() for s in sessions])
print(f"IC response (spikes/s): min {allR.min():.1f}, median {np.median(allR):.1f}, max {allR.max():.1f}")
print("Session-areas with >= 50 RS units:", {a: int((units_tab[f"{a}_RS"] >= 50).sum()) for a in il.AREAS},
      "(paper 12, 8, 11, 10, 12, 11)")
subject  V1_RS  LM_RS  RL_RS  AL_RS  PM_RS  AM_RS trials IC1/LC1/LC2/IC2 TRE1/TRE2 XRE1/XRE2
 619293     89      0    156    106    151    107        400/400/400/400   400/400   400/400
 619296    168      4     87     94     95     93        400/400/400/400   400/400   400/400
 620333    114    119    128     91    102    129        400/400/400/400   400/400   400/400
 620334     68    183    100      0     90     82        400/400/400/400   400/400   400/400
 625545    190    181    122    160    173    134        400/400/400/400   400/400   400/400
 625554    179      0    136    139    193    183        400/400/400/400   400/400   400/400
 625555    314      0    208    159    214    198        400/400/400/400   400/400   400/400
 630506     58    131    140    121    208      0        400/400/400/400   400/400   400/400
 631510    129    105     75     87     70     96        400/400/400/400   400/400   400/400
 631570    142     74      0    105    158    109        400/400/400/400   400/400   400/400
 633229    181    172    130    179    210    142        400/400/400/400   400/400   400/400
 637484    172    122    279      0     67     96        400/400/400/400   400/400   400/400
IC response (spikes/s): min 0.0, median 0.0, max 525.0
Session-areas with >= 50 RS units: {'V1': 12, 'LM': 8, 'RL': 11, 'AL': 10, 'PM': 12, 'AM': 11} (paper 12, 8, 11, 10, 12, 11)

The counts are those of the data check of decoding_shin_2025.ipynb.

Cross-validation seeds (fig3c, fig3d)¶

The tenfold partition is random. Reading (a) is repeated on every area with two more seeds (the seed of decoding_shin_2025.ipynb plus 100,000 and 200,000), the TRE decoder only. Reading (b) with the same seeds, on every area, is in decoding_seed1_shin_2025.ipynb and decoding_seed2_shin_2025.ipynb. The values of the first seed are those of the decoding notebook.

In [3]:
t0 = time.time()
jobs = []
for i, s in enumerate(sessions):
    u = s["units"]
    for j, a in enumerate(il.AREAS):
        m = ((u.area == a) & u.rs).values
        if m.sum() < 50:
            continue
        b = s["blocks"]["ICwcfg1"]
        jobs.append((s["subject"], a, b["R"][:, m], b["codes"], 1000 * i + j))

def tre(subject, a, R, codes, seed, rd, kernel="linear"):
    r = il.decode(R, codes, [il.IC1, il.LC1, il.LC2, il.IC2], [il.TRE1, il.TRE2], seed=seed, reading=rd, kernel=kernel)
    score, icv, lcv = il.inference_score(r["probe"])
    return {"subject": subject, "area": a, "seed": seed, "reading": rd, "kernel": kernel, "acc": r["acc"],
            "tre_ic": icv, "tre_lc": lcv, "inference": score, "unconverged": r["n_unconverged"]}

tasks = [(sub, a, R, c, sd + k * 100000, "a") for k in [1, 2] for (sub, a, R, c, sd) in jobs]
seeds = pd.DataFrame(Parallel(n_jobs=N_JOBS)(delayed(tre)(*t) for t in tasks))
seeds["seed_k"] = seeds.seed // 100000
print(f"{len(tasks)} decoders in {time.time() - t0:.0f} s; fits that did not converge: {int(seeds.unconverged.sum())}")
128 decoders in 560 s; fits that did not converge: 0
In [4]:
for (rd, k), x in seeds[seeds.area == "V1"].groupby(["reading", "seed_k"]):  # reading (a); (b) in decoding_seed1 and seed2
    print(f"V1, reading ({rd}), seed +{k}00000: P(TRE -> IC) {x.tre_ic.mean():.3f}, P(TRE -> LC) {x.tre_lc.mean():.3f}, "
          f"accuracy {x.acc.mean():.3f}; signed-rank IC vs LC p = {il.signrank_p(x.tre_ic, x.tre_lc):.4g}; "
          f"IC within 0.51 +/- 0.051: {abs(x.tre_ic.mean() - 0.51) <= 0.051}, LC within 0.22 +/- 0.022: {abs(x.tre_lc.mean() - 0.22) <= 0.022}")
V1, reading (a), seed +100000: P(TRE -> IC) 0.513, P(TRE -> LC) 0.198, accuracy 0.700; signed-rank IC vs LC p = 0.002441; IC within 0.51 +/- 0.051: True, LC within 0.22 +/- 0.022: True
V1, reading (a), seed +200000: P(TRE -> IC) 0.511, P(TRE -> LC) 0.199, accuracy 0.697; signed-rank IC vs LC p = 0.002441; IC within 0.51 +/- 0.051: True, LC within 0.22 +/- 0.022: True
In [5]:
rows = []
for k, x in seeds[seeds.reading == "a"].groupby("seed_k"):
    for a in il.AREAS:
        y = x[x.area == a]
        rows.append({"seed_k": k, "area": a, "sessions": len(y), "inference_mean": y.inference.mean(), "p": il.signrank_p(y.inference)})
t = pd.DataFrame(rows)
print(t.pivot(index="area", columns="seed_k", values=["inference_mean", "p"]).reindex(il.AREAS).round(4).to_string())
for k, y in t.groupby("seed_k"):
    sig = y[(y.p < 0.05) & (y.inference_mean > 0)].area.tolist()
    print(f"reading (a), seed +{k}00000: significant {sig} ({len(sig)} areas); three largest mean scores "
          f"{y.nlargest(3, 'inference_mean').area.tolist()}; three smallest p {y.nsmallest(3, 'p').area.tolist()}")
       inference_mean               p        
seed_k              1       2       1       2
area                                         
V1             0.3151  0.3119  0.0024  0.0024
LM             0.3031  0.3010  0.0156  0.0156
RL             0.0943  0.0922  0.0020  0.0029
AL             0.2572  0.2570  0.0020  0.0020
PM             0.0582  0.0563  0.0923  0.0923
AM             0.0620  0.0626  0.0068  0.0020
reading (a), seed +100000: significant ['V1', 'LM', 'RL', 'AL', 'AM'] (5 areas); three largest mean scores ['V1', 'LM', 'AL']; three smallest p ['RL', 'AL', 'V1']
reading (a), seed +200000: significant ['V1', 'LM', 'RL', 'AL', 'AM'] (5 areas); three largest mean scores ['V1', 'LM', 'AL']; three smallest p ['AL', 'AM', 'V1']

Result: in V1, under reading (a), both TRE values are within tolerance with each further seed (IC 0.513 and 0.511, LC 0.198 and 0.199), with signed-rank p = 0.0024 for both. Each further seed gives the same five significant areas as the decoding notebook (V1, LM, RL, AL, AM; PM p = 0.092) and the same three largest mean scores (V1, LM, AL). Under reading (a), the fig3c and fig3d results do not depend on the seed; reading (b) is in the two seed notebooks.

Bound of the C search and coding (decoder reading b)¶

MATLAB's OptimizeHyperparameters for fitcecoc searches the coding (one-vs-one or one-vs-all) and the SVM's box constraint and kernel scale over 10^-3 to 10^3. For a linear kernel, the kernel scale s divides the inputs, which gives the same classifier as C / s^2, so a search over C covers both. Reading (b) searches C over 10^-3 to 1 with one-vs-one coding. Here, on the training trials of the first outer fold of each session, in V1 and in PM and AM (the two areas whose significance in fig3d depends on the decoder), the inner threefold accuracy is computed for C from 10^-3 to 10^3 (one-vs-one) and for C up to 1 with one-vs-all coding, to show whether the bound and the coding can change the choice.

In [6]:
t0 = time.time()
path_jobs = [j for j in jobs if j[1] in ("V1", "PM", "AM")]
cp = Parallel(n_jobs=N_JOBS)(delayed(il.c_path)(R, c, [il.IC1, il.LC1, il.LC2, il.IC2], sd) for (sub, a, R, c, sd) in path_jobs)
cp = pd.concat([d.assign(subject=sub, area=a) for d, (sub, a, *_ ) in zip(cp, path_jobs)])
print(f"computed in {time.time() - t0:.0f} s")
for a in ["V1", "PM", "AM"]:
    x = cp[cp.area == a]
    print(f"\n{a}, inner threefold accuracy by coding and C:")
    print(x.pivot_table(index="subject", columns=["coding", "C"], values="inner_cv_acc").round(3).to_string())
    best_le1 = x[x.C <= 1].groupby("subject").inner_cv_acc.max()
    best_gt1 = x[x.C > 1].groupby("subject").inner_cv_acc.max()
    best_ovo = x[(x.C <= 1) & (x.coding == "ovo")].groupby("subject").inner_cv_acc.max()
    best_ovr = x[(x.C <= 1) & (x.coding == "ovr")].groupby("subject").inner_cv_acc.max()
    d = best_gt1 - best_le1
    print(f"{a}: best with C > 1 minus best with C <= 1 (both codings), per session:", d.round(3).to_dict())
    print(f"{a}: sessions where C > 1 is better: {int((d > 0).sum())} of {len(d)}; range {d.min():.3f} to {d.max():.3f}")
    e = best_ovr - best_ovo
    print(f"{a}: best one-vs-all minus best one-vs-one, C <= 1: range {e.min():.3f} to {e.max():.3f}")
    print(f"{a}: one-vs-one fits that did not converge, by C (of {3 * x.subject.nunique()} per C):",
          x[x.coding == "ovo"].groupby("C").unconverged.sum().to_dict())
computed in 792 s

V1, inner threefold accuracy by coding and C:
coding       ovo                                                            ovr                           
C       0.001    0.010    0.100    1.000    10.000   100.000  1000.000 0.001    0.010    0.100    1.000   
subject                                                                                                   
619293     0.253    0.267    0.253    0.263    0.256    0.271    0.267    0.267    0.276    0.254    0.253
619296     0.908    0.915    0.895    0.898    0.898    0.898    0.898    0.908    0.906    0.897    0.890
620333     0.652    0.696    0.688    0.673    0.662    0.658    0.605    0.638    0.696    0.689    0.687
620334     0.588    0.685    0.685    0.676    0.671    0.669    0.620    0.672    0.690    0.710    0.699
625545     0.803    0.811    0.776    0.751    0.751    0.751    0.751    0.798    0.815    0.799    0.763
625554     0.773    0.806    0.775    0.767    0.760    0.760    0.760    0.772    0.804    0.783    0.765
625555     0.509    0.510    0.497    0.503    0.499    0.499    0.499    0.497    0.498    0.481    0.483
630506     0.635    0.731    0.746    0.730    0.726    0.718    0.697    0.654    0.708    0.730    0.734
631510     0.823    0.851    0.833    0.811    0.807    0.804    0.797    0.822    0.849    0.853    0.835
631570     0.790    0.812    0.798    0.767    0.768    0.768    0.768    0.810    0.827    0.799    0.774
633229     0.860    0.873    0.866    0.845    0.845    0.845    0.845    0.881    0.881    0.862    0.842
637484     0.544    0.563    0.565    0.557    0.556    0.558    0.528    0.553    0.564    0.559    0.553
V1: best with C > 1 minus best with C <= 1 (both codings), per session: {'619293': -0.005, '619296': -0.017, '620333': -0.033, '620334': -0.039, '625545': -0.064, '625554': -0.046, '625555': -0.011, '630506': -0.02, '631510': -0.047, '631570': -0.059, '633229': -0.036, '637484': -0.006}
V1: sessions where C > 1 is better: 0 of 12; range -0.064 to -0.005
V1: best one-vs-all minus best one-vs-one, C <= 1: range -0.012 to 0.024
V1: one-vs-one fits that did not converge, by C (of 36 per C): {0.001: 0, 0.01: 0, 0.1: 0, 1.0: 0, 10.0: 12, 100.0: 17, 1000.0: 17}

PM, inner threefold accuracy by coding and C:
coding       ovo                                                            ovr                           
C       0.001    0.010    0.100    1.000    10.000   100.000  1000.000 0.001    0.010    0.100    1.000   
subject                                                                                                   
619293     0.266    0.267    0.271    0.259    0.267    0.270    0.259    0.260    0.248    0.248    0.244
619296     0.544    0.558    0.558    0.567    0.564    0.555    0.468    0.520    0.557    0.555    0.561
620333     0.331    0.374    0.373    0.362    0.360    0.354    0.294    0.349    0.359    0.373    0.370
620334     0.606    0.614    0.604    0.600    0.588    0.591    0.534    0.598    0.610    0.609    0.619
625545     0.362    0.407    0.397    0.378    0.358    0.334    0.328    0.404    0.424    0.408    0.404
625554     0.712    0.696    0.669    0.646    0.649    0.649    0.649    0.691    0.706    0.681    0.655
625555     0.570    0.599    0.556    0.510    0.486    0.486    0.486    0.572    0.576    0.551    0.546
630506     0.519    0.567    0.535    0.521    0.512    0.517    0.498    0.567    0.556    0.556    0.539
631510     0.356    0.440    0.432    0.439    0.429    0.429    0.362    0.426    0.424    0.431    0.438
631570     0.644    0.701    0.703    0.686    0.695    0.692    0.670    0.657    0.697    0.697    0.698
633229     0.521    0.549    0.511    0.499    0.493    0.485    0.480    0.528    0.548    0.518    0.497
637484     0.346    0.379    0.392    0.394    0.394    0.394    0.319    0.356    0.369    0.364    0.358
PM: best with C > 1 minus best with C <= 1 (both codings), per session: {'619293': -0.001, '619296': -0.003, '620333': -0.014, '620334': -0.028, '625545': -0.066, '625554': -0.063, '625555': -0.113, '630506': -0.051, '631510': -0.01, '631570': -0.008, '633229': -0.056, '637484': 0.001}
PM: sessions where C > 1 is better: 1 of 12; range -0.113 to 0.001
PM: best one-vs-all minus best one-vs-one, C <= 1: range -0.024 to 0.017
PM: one-vs-one fits that did not converge, by C (of 36 per C): {0.001: 0, 0.01: 0, 0.1: 0, 1.0: 0, 10.0: 28, 100.0: 30, 1000.0: 30}

AM, inner threefold accuracy by coding and C:
coding       ovo                                                            ovr                           
C       0.001    0.010    0.100    1.000    10.000   100.000  1000.000 0.001    0.010    0.100    1.000   
subject                                                                                                   
619293     0.258    0.249    0.254    0.242    0.255    0.244    0.258    0.267    0.262    0.285    0.274
619296     0.408    0.475    0.455    0.451    0.456    0.446    0.331    0.466    0.467    0.457    0.460
620333     0.351    0.405    0.385    0.383    0.385    0.383    0.316    0.406    0.409    0.421    0.413
620334     0.485    0.525    0.508    0.497    0.490    0.490    0.455    0.499    0.494    0.497    0.497
625545     0.463    0.525    0.510    0.489    0.483    0.481    0.421    0.497    0.500    0.504    0.496
625554     0.665    0.661    0.632    0.610    0.612    0.612    0.612    0.642    0.663    0.622    0.610
625555     0.460    0.475    0.459    0.441    0.430    0.426    0.428    0.468    0.479    0.464    0.458
631510     0.448    0.496    0.481    0.469    0.465    0.465    0.408    0.506    0.500    0.494    0.495
631570     0.505    0.536    0.522    0.496    0.494    0.493    0.447    0.497    0.511    0.526    0.512
633229     0.571    0.604    0.558    0.532    0.516    0.513    0.487    0.563    0.583    0.571    0.572
637484     0.284    0.298    0.303    0.305    0.300    0.300    0.278    0.297    0.289    0.299    0.292
AM: best with C > 1 minus best with C <= 1 (both codings), per session: {'619293': -0.026, '619296': -0.019, '620333': -0.036, '620334': -0.035, '625545': -0.042, '625554': -0.052, '625555': -0.049, '631510': -0.04, '631570': -0.042, '633229': -0.088, '637484': -0.005}
AM: sessions where C > 1 is better: 0 of 11; range -0.088 to -0.005
AM: best one-vs-all minus best one-vs-one, C <= 1: range -0.026 to 0.027
AM: one-vs-one fits that did not converge, by C (of 33 per C): {0.001: 0, 0.01: 0, 0.1: 0, 1.0: 0, 10.0: 30, 100.0: 30, 1000.0: 30}

Result: the best inner accuracy with C above 1 is lower than the best with C up to 1 in every V1 and AM session (by 0.005 to 0.064 in V1 and 0.005 to 0.088 in AM) and in 11 of 12 PM sessions (by 0.001 to 0.113); in PM session 637484 it is higher by 0.001. Most fits with C above 1 stop at max_iter (12 to 17 of 36 per C value in V1, 28 to 30 of 36 in PM, 30 of 33 in AM), so these accuracies come partly from unconverged fits, and the check covers the first outer fold only. On these fits, extending the grid above C = 1 would change the chosen C in at most one session of the three areas, by 0.001 in inner accuracy. The one-vs-all coding changes the best inner accuracy by -0.012 to 0.024 in V1, -0.024 to 0.017 in PM and -0.026 to 0.027 in AM, so it could change the choice in some sessions; its effect on the inference scores was not computed.

supp3-decoder-types¶

Decoders of other kinds (quadratic and radial basis function SVMs, a fully connected neural network) give results highly consistent with the linear SVM.

The quadratic and RBF SVMs replace the linear kernel in decoder reading (a) (C = 1), with the first extra seed of the section above, on V1, LM and AL; the paper gives no kernel parameters, so scikit-learn's data-scaled gamma="scale" is assumed (helper, decode). Each is compared with the linear SVM on the same seed: sign of the mean TRE inference score and significance of the signed-rank test against 0 (p < 0.05). The network of the Methods (two hidden layers of 100 units, softmax output, Adam, cross-entropy, 50,000 epochs, tenfold) is timed on one fold of one session, and its full cost estimated.

In [7]:
t0 = time.time()
tasks = [(sub, a, R, c, sd + 100000, "a", k) for k in ["poly", "rbf"] for (sub, a, R, c, sd) in jobs if a in ("V1", "LM", "AL")]
kern = pd.DataFrame(Parallel(n_jobs=N_JOBS)(delayed(tre)(*t) for t in tasks))
kern = pd.concat([kern, seeds[(seeds.reading == "a") & (seeds.seed_k == 1) & seeds.area.isin(["V1", "LM", "AL"])]])
print(f"{len(tasks)} decoders in {time.time() - t0:.0f} s")
rows = []
for (a, k), x in kern.groupby(["area", "kernel"]):
    rows.append({"area": a, "kernel": k, "sessions": len(x), "accuracy": x.acc.mean(), "IC": x.tre_ic.mean(), "LC": x.tre_lc.mean(),
                 "inference_mean": x.inference.mean(), "p": il.signrank_p(x.inference)})
kt = pd.DataFrame(rows)
kt["sig_positive"] = (kt.p < 0.05) & (kt.inference_mean > 0)
print(kt.round(4).to_string(index=False))
lin = kt[kt.kernel == "linear"].set_index("area")
for k in ["poly", "rbf"]:
    y = kt[kt.kernel == k].set_index("area")
    agree = {a: bool((np.sign(y.loc[a, "inference_mean"]) == np.sign(lin.loc[a, "inference_mean"])) and
                     (y.loc[a, "sig_positive"] == lin.loc[a, "sig_positive"])) for a in ["V1", "LM", "AL"]}
    print(f"{k}: same sign and significance as linear in", agree)
60 decoders in 42 s
area kernel  sessions  accuracy     IC     LC  inference_mean      p  sig_positive
  AL linear        10    0.6491 0.5108 0.2537          0.2572 0.0020          True
  AL   poly        10    0.6879 0.5302 0.2610          0.2692 0.0039          True
  AL    rbf        10    0.6762 0.5250 0.2627          0.2623 0.0020          True
  LM linear         8    0.7074 0.5379 0.2348          0.3031 0.0156          True
  LM   poly         8    0.7451 0.5626 0.2279          0.3347 0.0156          True
  LM    rbf         8    0.7315 0.5446 0.2386          0.3060 0.0156          True
  V1 linear        12    0.6998 0.5132 0.1981          0.3151 0.0024          True
  V1   poly        12    0.7301 0.5295 0.2106          0.3190 0.0024          True
  V1    rbf        12    0.7201 0.4903 0.2179          0.2724 0.0015          True
poly: same sign and significance as linear in {'V1': True, 'LM': True, 'AL': True}
rbf: same sign and significance as linear in {'V1': True, 'LM': True, 'AL': True}
In [8]:
from sklearn.neural_network import MLPClassifier
s = sessions[0]
b = s["blocks"]["ICwcfg1"]
m = ((s["units"].area == "V1") & s["units"].rs).values
sel = np.isin(b["codes"], [il.IC1, il.LC1, il.LC2, il.IC2])
X = stats.zscore(b["R"][sel][:, m], axis=0)
X[~np.isfinite(X)] = 0
y = b["codes"][sel]
n_ep = 300
net = MLPClassifier(hidden_layer_sizes=(100, 100), solver="adam", batch_size=len(y), max_iter=n_ep, tol=0,
                    n_iter_no_change=10 ** 9, alpha=0.0, random_state=0)
t0 = time.time()
net.fit(X[: int(0.9 * len(y))], y[: int(0.9 * len(y))])
per_epoch = (time.time() - t0) / n_ep
n_fits = 10 * len([j for j in jobs if j[1] in ("V1", "LM", "AL")])
print(f"{per_epoch * 1000:.2f} ms per full-batch epoch on {m.sum()} units; one fit of 50,000 epochs ~ {per_epoch * 5e4 / 60:.1f} min; "
      f"{n_fits} fits (10 folds x {n_fits // 10} session-areas of V1, LM, AL) ~ {per_epoch * 5e4 * n_fits / 3600:.1f} h on one core")
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/neural_network/_multilayer_perceptron.py:695: UserWarning: Got `batch_size` less than 1 or larger than sample size. It is going to be clipped
  warnings.warn(
4.38 ms per full-batch epoch on 89 units; one fit of 50,000 epochs ~ 3.6 min; 300 fits (10 folds x 30 session-areas of V1, LM, AL) ~ 18.2 h on one core
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/neural_network/_multilayer_perceptron.py:785: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (300) reached and the optimization hasn't converged yet.
  warnings.warn(

Result: blocked (too_large) for the network; the two SVMs agree with the linear SVM. With the quadratic and RBF kernels, the TRE inference score is positive and significant in V1, LM and AL, as with the linear SVM on the same seed (p from 0.0015 to 0.0156 for all three kernels). The network of the Methods (50,000 full-batch epochs, tenfold) takes the time per epoch printed above, 2.5 to 4.4 ms in the runs so far, which is 2 to 4 minutes per fit and about 11 to 18 hours on one core for the 300 fits of these three areas, far beyond the notebook's limit, so the claim, which needs each decoder type, is not decided.

In [9]:
print(f"Notebook ran in {(time.time() - started) / 60:.1f} min")
Notebook ran in 25.5 min

Summary¶

Check or claim Result
Seeds for fig3c and fig3d, reading (a) fig3c within tolerance and fig3d five significant areas with each seed; reading (b) in the seed notebooks
Bound of the C search and coding (V1, PM, AM) C above 1 better in 1 of 35 sessions, by 0.001, from mostly unconverged fits; one-vs-all changes the best inner accuracy by -0.026 to 0.027
supp3-decoder-types blocked (too_large) for the network; quadratic and RBF SVMs agree with the linear SVM in V1, LM and AL

How the analysis differs from the paper's: the quadratic and RBF kernels use scikit-learn's data-scaled gamma="scale", as the paper gives no kernel parameters; the network is timed with scikit-learn's MLPClassifier, where the paper used PyTorch.