Spaces:
Running on Zero
Running on Zero
RootScope v4: manuscript model (LightGBM per round, 3 seeds, layer + soft-neighbor context, radial prior)
47c4bf8 verified | """Soft radial-order constraint (fix for error family B: one-layer boundary shifts). | |
| POST-HOC. Takes the model's probability matrix and the cell centroids; it never looks at | |
| the 768 DINOv2 columns, so it can be re-applied unchanged after the re-fine-tune. | |
| Method | |
| 1. Per section, a monotone radial coordinate r_rel in [0,1]: r_rel = r / r_surface(theta), | |
| centre = mean centroid of the cells PREDICTED vascular-core, r_surface = furthest cell | |
| within +-15 deg of the cell's own angle (all cells of the section, labels not needed). | |
| 2. Tissues get a radial RANK: epidermis/root_cap 0 < exodermis 1 < cortex 2 < endodermis 3 | |
| < pericycle 4 < {stele,xylem,phloem} 5. The trio shares ONE rank - there is no reliable | |
| radial order inside it, so the constraint cannot touch stele/xylem/phloem against each other. | |
| 3. Per section, an ordered-logit model of rank against -r_rel is fitted to the model's OWN | |
| posterior (soft targets, no labels, nothing from the ground truth). Two passes: the second | |
| fit uses the once-corrected posterior. | |
| 4. P_final ~ P_model * P_order ** lambda, renormalised. | |
| Hyperparameters (s = 0.3, lambda = 2.0, 2 passes) were chosen on the 100 out-of-fold training | |
| sections only; the OOF optimum is a broad plateau (s 0.3-1.2, 2-4 passes, lambda 2). | |
| Usage | |
| from apply_radial_order import correct | |
| P_final = correct(P, source_file, cell_id, classes) # P: (n_cells, n_classes) | |
| """ | |
| import numpy as np | |
| from scipy.optimize import minimize | |
| CLASSES_DEFAULT = ['cortex','endodermis','epidermis','exodermis','pericycle','phloem','root_cap','stele','xylem'] | |
| RANK = {'epidermis':0,'root_cap':0,'exodermis':1,'cortex':2,'endodermis':3,'pericycle':4, | |
| 'stele':5,'xylem':5,'phloem':5} | |
| NRANK = 6 | |
| CORE = ('endodermis','pericycle','stele','xylem','phloem') | |
| S_SCALE, LAMBDA, PASSES, FLOOR = 0.3, 2.0, 2, 1e-3 | |
| def radial_coordinate(cy, cx, cy_all, cx_all, pred_names, win_deg=15.0): | |
| """r_rel for the cells (cy,cx); the surface envelope uses every cell of the section | |
| (cy_all,cx_all), which may be a superset (unlabelled cells included).""" | |
| core = np.isin(pred_names, CORE) | |
| if core.sum() >= 5: | |
| y0, x0 = cy[core].mean(), cx[core].mean() | |
| else: | |
| y0, x0 = cy_all.mean(), cx_all.mean() | |
| r_all = np.hypot(cy_all - y0, cx_all - x0); th_all = np.arctan2(cy_all - y0, cx_all - x0) | |
| r = np.hypot(cy - y0, cx - x0); th = np.arctan2(cy - y0, cx - x0) | |
| da = np.abs((th[:, None] - th_all[None, :] + np.pi) % (2 * np.pi) - np.pi) | |
| rs = np.maximum(np.where(da < np.deg2rad(win_deg), r_all[None, :], -1.0).max(1), np.maximum(r, 1.0)) | |
| return r / rs | |
| def _nll(par, v, Q, s): | |
| t = np.concatenate([[par[0]], par[0] + np.cumsum(np.logaddexp(0.0, par[1:]))]) | |
| z = np.clip((t[None, :] - v[:, None]) / s, -50, 50) | |
| cum = np.concatenate([np.zeros((len(v), 1)), 1 / (1 + np.exp(-z)), np.ones((len(v), 1))], 1) | |
| return -(Q * np.log(np.clip(np.diff(cum, axis=1), 1e-9, 1.0))).sum() | |
| def _order_prior(v, Q, s): | |
| cw = np.clip(np.cumsum(Q.sum(0) + 1e-9)[:-1] / (Q.sum() + 1e-9), 0.01, 0.99) | |
| t0 = np.maximum.accumulate(np.quantile(v, cw)) | |
| par0 = np.concatenate([[t0[0]], np.log(np.expm1(np.maximum(np.diff(t0), 1e-3)))]) | |
| par = minimize(_nll, par0, args=(v, Q, s), method='L-BFGS-B', options=dict(maxiter=500)).x | |
| t = np.concatenate([[par[0]], par[0] + np.cumsum(np.logaddexp(0.0, par[1:]))]) | |
| z = np.clip((t[None, :] - v[:, None]) / s, -50, 50) | |
| cum = np.concatenate([np.zeros((len(v), 1)), 1 / (1 + np.exp(-z)), np.ones((len(v), 1))], 1) | |
| return np.clip(np.diff(cum, axis=1), 1e-9, 1.0), t | |
| def correct_section(P, cy, cx, classes, cy_all=None, cx_all=None, | |
| s=S_SCALE, lam=LAMBDA, passes=PASSES): | |
| cy_all = cy if cy_all is None else cy_all; cx_all = cx if cx_all is None else cx_all | |
| pred = np.asarray(classes)[P.argmax(1)] | |
| v = -radial_coordinate(cy, cx, cy_all, cx_all, pred) | |
| Pc = P | |
| for k in range(passes): | |
| Q = np.zeros((len(P), NRANK)) | |
| for j, c in enumerate(classes): | |
| Q[:, RANK[c]] += Pc[:, j] | |
| Q /= np.maximum(Q.sum(1, keepdims=True), 1e-12) | |
| Po, _ = _order_prior(v, Q, s) | |
| Po = Po + FLOOR; Po /= Po.sum(1, keepdims=True) | |
| pri = np.stack([Po[:, RANK[c]] for c in classes], 1) | |
| F = P * pri ** lam; Pc = F / F.sum(1, keepdims=True) | |
| return Pc | |
| def correct(P, source_file, cell_id, classes=None, geom=None, **kw): | |
| """P (n,k) posterior; source_file (n,) section of each row. geom: optional dict | |
| section -> (cy_all, cx_all, cy, cx) if the centroids are not in geom.npz.""" | |
| classes = classes or CLASSES_DEFAULT | |
| if geom is None: | |
| raise ValueError( | |
| "correct() needs cell centroids. Pass geom={source_file: (cy_all, cx_all, " | |
| "{cell_id: (cy, cx)})}. The research build read these from a cached " | |
| "geom.npz that is not part of this package.") | |
| d = None # unreachable; kept so the original block below still parses | |
| gsf = d['source_file'].astype(str); gid = d['cell_id'].astype(np.int64) | |
| geom = {} | |
| for s_ in np.unique(source_file): | |
| m = gsf == s_ | |
| geom[s_] = (d['cy'][m], d['cx'][m], dict(zip(gid[m], zip(d['cy'][m], d['cx'][m])))) | |
| out = P.copy() | |
| for s_ in np.unique(source_file): | |
| m = source_file == s_ | |
| cy_all, cx_all, look = geom[s_] | |
| yx = np.array([look[int(c)] for c in np.asarray(cell_id)[m]]) | |
| out[m] = correct_section(P[m], yx[:, 0], yx[:, 1], classes, cy_all, cx_all, **kw) | |
| return out | |