Optimal coloring (Yemini et al. 2021): other sparsities for 3 reporters

Paper. Yemini et al. (2021), Cell 184:272-288, doi:10.1016/j.cell.2020.12.012.

Data. The neighbors and reporters of optimal_coloring.ipynb: the head atlas fit on the ten heads of DANDI:000715, version 0.241009.1514, and the Brain Atlas reporters of data.mat of the paper's Optimal_Coloring code (github.com/Eviatar/Optimal_Coloring at 990c597), an outside input.

The sparsity that gives the 3-reporter solution of fig7de-optimal-coloring-violations is not given. optimal_coloring.ipynb runs the criterion's 7,500 runs at 1500, the sparsity of its scan whose runs most often keep 3 reporters. This notebook runs the scan's other sparsities whose runs keep 3 reporters, apart so that each notebook stays within its time limit: 2000 and 3000, where most runs keep 3 reporters in the scan, with the criterion's 7,500 runs, and 1000, where about half do, with 1,000 runs. The conclusion is in optimal_coloring.ipynb.

import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import neuropal as npl

pd.set_option("display.width", 200)
started = time.time()
worms = npl.load_all()
oc = npl.load_optimal_coloring()
print(f"loaded {len(worms)} files from dandiset {npl.DANDISET_ID} version {npl.DANDISET_VERSION} and "
      f"Optimal_Coloring data.mat in {time.time() - started:.0f} s")

loaded 10 files from dandiset 000715 version 0.241009.1514 and Optimal_Coloring data.mat in 18 s

Data check

feats = [npl.worm_features(w) for w in worms]
print("annotated neurons per worm:", {w["subject"]: len(f[0]) for w, f in zip(worms, feats)})
names, mu, sigma, atlas_info = npl.fit_atlas(feats, seed=0)
print(f"head atlas: {len(names)} neurons, worms per neuron {atlas_info['counts'].min()} to {atlas_info['counts'].max()}, "
      f"positions {mu[:, :3].min():.1f} to {mu[:, :3].max():.1f} um")
BA = oc["BrainAtlas"]
print(f"\ndata.mat BrainAtlas table: {BA.shape[0]} neurons x {BA.shape[1] - 3} reporter columns (after 3 descriptive columns "
      f"{list(BA.columns[:3])}); values {sorted(np.unique(BA.iloc[:, 3:].to_numpy()).tolist())}")
print(f"data.mat neurons to color: {len(oc['neurons'])}; all in the table: {set(oc['neurons']) <= set(BA.Neuron.astype(str))}")
R_ref, rep_ref = npl.reporter_matrix(BA, oc["neurons"])
print(f"reporters with a distinct pattern over these {len(oc['neurons'])} neurons: {R_ref.shape[1]} "
      f"(Methods S2 section 4.1: 1192 reporters, unique = 654)")
A_ref = np.asarray(oc["neuron_prob_distances"]).astype(int)
print(f"data.mat adjacency: {A_ref.shape}, values {np.unique(A_ref).tolist()}, symmetric {bool((A_ref == A_ref.T).all())}, "
      f"diagonal entries set {int(np.diag(A_ref).sum())}, adjacent ordered pairs {int(A_ref.sum())}")

annotated neurons per worm: {'11_YAaLR': 193, '20_YAaLR': 193, '38_YAaDV': 190, '55_YAaDV': 192, '56_YAaDV': 191, '62_YAaLR': 186, '64_YAaDV': 189, '7_YAaLR': 188, '70_YAaLR': 193, '76_YAaDV': 188}
head atlas: 198 neurons, worms per neuron 1 to 10, positions 5.4 to 187.1 um

data.mat BrainAtlas table: 302 neurons x 1189 reporter columns (after 3 descriptive columns ['NeuronClass', 'Neuron', 'NeuronType']); values [0.0, 1.0]
data.mat neurons to color: 235; all in the table: True
reporters with a distinct pattern over these 235 neurons: 654 (Methods S2 section 4.1: 1192 reporters, unique = 654)
data.mat adjacency: (235, 235), values [0, 1], symmetric False, diagonal entries set 0, adjacent ordered pairs 5655

# The neighbors of ms2-neuropal-color-violations: atlas neurons that data.mat colors, adjacency from
# computeProbAdjacency on the head atlas.
shared = [n for n in names if n in oc["neurons"]]
ia = [names.index(n) for n in shared]
A_head = npl.prob_adjacency(mu[ia, :3], sigma[ia][:, :3, :3])
np.fill_diagonal(A_head, 0)
R_head, rep_head = npl.reporter_matrix(BA, shared)
print(f"head atlas neurons in data.mat: {len(shared)} of {len(names)}; left out: {[n for n in names if n not in oc['neurons']]}")
print(f"head adjacency: {int(A_head.sum())} adjacent ordered pairs, symmetric {bool((A_head == A_head.T).all())}, "
      f"neighbors per neuron {A_head.sum(1).min()} to {A_head.sum(1).max()} (median {int(np.median(A_head.sum(1)))})")
print(f"reporters with a distinct pattern over the {len(shared)} head neurons: {R_head.shape[1]}")
j = [oc["neurons"].index(n) for n in shared]
print(f"agreement with data.mat's adjacency on the same neurons: {100 * (A_head == A_ref[np.ix_(j, j)]).mean():.1f}% of entries")

head atlas neurons in data.mat: 192 of 198; left out: ['AMSOL', 'AMSOR', 'DA2', 'DB3', 'HMC', 'VD3']
head adjacency: 4065 adjacent ordered pairs, symmetric False, neighbors per neuron 0 to 70 (median 16)
reporters with a distinct pattern over the 192 head neurons: 543
agreement with data.mat's adjacency on the same neurons: 96.1% of entries

fig7de-optimal-coloring-violations: other sparsities

The runs, best-run selection and seeds are those of optimal_coloring.ipynb.

K, MARGIN, ITERS, BATCH = 3, 1 / 3, 300, 250

def run_many(A, R, sparsity, runs, rng):
    out = []
    for b in range(0, runs, BATCH):
        X = npl.optimal_color_solver(A, R, K, MARGIN, sparsity, ITERS, min(BATCH, runs - b), rng)
        out.append(npl.solution_violations(X, R, A, MARGIN))
    t = pd.DataFrame(np.concatenate(out), columns=["reporters", "share", "shortfall"])
    t["reporters"] = t.reporters.astype(int)
    t[["share", "shortfall"]] *= 100
    return t

def best(t, label, reporters=None):
    s = t if reporters is None else t[t.reporters == reporters]
    if len(s) == 0:
        return dict(selection=label, runs=0)
    r = s.loc[s.shortfall.idxmin()]
    return dict(selection=label, runs=len(s), reporters_of_best=int(r.reporters), shortfall_pct=r.shortfall,
                methods_share_pct=r.share, lowest_share_pct_any_run=s.share.min())

t0 = time.time()
rows = []
for k, (s, n) in enumerate([(1000, 1_000), (2000, 7_500), (3000, 7_500)]):
    t = run_many(A_head, R_head, s, n, np.random.default_rng(7146 + k))
    v3 = t[t.reporters == 3].shortfall
    rows.append(dict(sparsity=s, **best(t, "best run with 3 reporters", 3), share_of_runs_with_3_reporters=(t.reporters == 3).mean(),
                     median_reporters=int(t.reporters.median()), p5_3_reporter_pct=v3.quantile(0.05), runs_below_9_45=int((v3 < 9.45).sum()),
                     distinct_values_up_to_10_pct=np.unique(v3[v3 <= 10].round(3)).tolist()))
print(f"runs in {time.time() - t0:.0f} s")
print(pd.DataFrame(rows).round(3).to_string(index=False))
print("\nWithin 9.45 to 11.55%:", {r["sparsity"]: bool(9.45 <= r["shortfall_pct"] <= 11.55) for r in rows})
print(f"notebook time {(time.time() - started) / 60:.1f} min")

runs in 1028 s
 sparsity                 selection  runs  reporters_of_best  shortfall_pct  methods_share_pct  lowest_share_pct_any_run  share_of_runs_with_3_reporters  median_reporters  p5_3_reporter_pct  runs_below_9_45 distinct_values_up_to_10_pct
     1000 best run with 3 reporters   561                  3          7.831             23.493                    23.493                           0.561                 3              9.897                2         [7.831, 8.38, 9.897]
     2000 best run with 3 reporters  7020                  3          8.380             25.141                    25.141                           0.936                 3              9.897                1  [8.38, 9.537, 9.717, 9.897]
     3000 best run with 3 reporters  6253                  3          9.897             29.692                    29.692                           0.834                 3             17.524                0                      [9.897]

Within 9.45 to 11.55%: {1000: False, 2000: False, 3000: True}
notebook time 17.5 min