You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Internal utilities for reproducible, per-batch cell subsampling.
27
-
`max_cells` caps pert cells only; ctrl cells are added on top.
28
-
- Both batch functions accept DataFrame `A` with named perturbation columns.
29
20
-`n_batches` parameter for `fit_gcate_batch` and `gcate_lfc_batch`:
30
21
specifies total number of batches instead of per-batch count; overrides
31
-
`batch_size` when set. Batches are sized evenly with `numpy.array_split`
32
-
so the last batch is never more than 1 perturbation smaller than the
33
-
others, avoiding an unstable single-pert tail batch.
22
+
`batch_size` when set.
34
23
-`estimate_r(max_cells=N, random_state=0)`: new parameter that
35
24
automatically subsamples to at most `N` cells before running JIC
36
-
selection, prioritising control cells (all-zero rows of `A`). Filling
37
-
the ctrl budget first and then drawing remaining slots from treated cells
38
-
is equivalent to full-data estimation because confounding structure is
39
-
concentrated in the baseline transcriptome.
25
+
selection, prioritising control cells.
40
26
41
-
### Fixed
42
-
43
-
-**Numba TBB fork warning**: Set `NUMBA_THREADING_LAYER_PRIORITY` to prefer OpenMP over TBB in `__init__.py`, eliminating `Numba: Attempted to fork from a non-main thread, the TBB library may be in an invalid state in the child process` warnings when Joblib forks after Numba parallel execution.
44
-
- Added `llvm-openmp` to conda dependencies for fork-safe OpenMP threading.
45
-
-**Fast-path threshold** (`gcate_glm.py`): Raised the effective design-dimension ceiling and added a throughput heuristic; GCATE runs with larger `r` or wide batch designs now correctly use the crispyx path.
46
-
-**Backend toggle** (`gcate_glm.py`): Re-added `_USE_FAST_BACKEND` module flag and `_backend_override()` context manager; setting `_USE_FAST_BACKEND = False` now reliably forces statsmodels through the entire GCATE + LFC call graph.
47
-
-**Weighted dispersion** (`nb_glm_fast.py`): Dispersion averaging is now cell-count-weighted; low-coverage perturbations contribute proportionally less to the pooled estimate.
48
-
-**Control-cell residuals** (`nb_glm_fast.py`): Control-cell deviance residuals and fitted values are now initialised from the global covariate model, eliminating a last-perturbation overwrite bug.
49
-
-**Module-qualified imports** (`gcate_opt.py`, `gcate.py`, `DR_estimation.py`): Replaced import-time name bindings with module-qualified references so backend toggles propagate correctly at call time.
50
-
-**`estimate_r` bare name** (`gcate.py`): Fixed a `NameError` when calling `estimate_r` that was caused by a bare `fit_glm_auto` reference.
51
-
-**crispyx availability check** (`gcate_glm.py`): Added a one-time import check; users without crispyx now get a transparent fallback to statsmodels instead of a traceback.
52
-
-**Dead code removal** (`nb_glm_fast.py`): Removed an unreachable branch in `_compute_poisson_deviance_residuals`.
53
-
54
-
### Added
55
-
56
-
-**Fast GLM backend via crispyx** (`nb_glm_fast.py`)
57
-
-`fit_glm_fast()`: Batch NB-GLM fitting using crispyx's `NBGLMBatchFitter`, replacing per-gene statsmodels IRLS with vectorized batch IRLS across all genes simultaneously.
-`fit_glm_ondisk()`: On-disk streaming GLM fitting for large h5ad files.
32
+
- Per-perturbation fitting (`_fit_glm_fast_per_perturbation`): for
33
+
multi-treatment data, fits binary (ctrl vs. treatment_k) models
34
+
independently, then assembles the full coefficient matrix.
35
+
-`fit_glm_auto()`: Routes to `fit_glm_fast()` when crispyx is available and
36
+
the effective design dimension is small; falls back to statsmodels otherwise.
37
+
-`estimate_disp_auto()`: Routes to `estimate_disp_fast()` for large gene
38
+
counts; falls back to statsmodels otherwise.
60
39
61
-
-**Per-perturbation fitting strategy** (`_fit_glm_fast_per_perturbation()` in `nb_glm_fast.py`)
62
-
- For multi-treatment data, loops over perturbation columns and fits binary (control vs. treatment_k) models, avoiding the memory and runtime cost of wide joint design matrices.
63
-
- Covariate coefficients are averaged across perturbation-specific models; treatment coefficients are assembled into the full coefficient matrix.
64
-
- Imputation is performed per-perturbation with correct counterfactual construction.
40
+
### Fixed
65
41
66
-
-**Auto-dispatch wrappers** (`gcate_glm.py`)
67
-
-`fit_glm_auto()`: Routes to `fit_glm_fast()` when the effective design dimension is small and crispyx is available; falls back to statsmodels otherwise. Includes convergence checking with automatic fallback.
68
-
-`estimate_disp_auto()`: Routes to `estimate_disp_fast()` when the number of genes is large enough for vectorized estimation to be faster.
42
+
-**Numba TBB fork warning**: Set `NUMBA_THREADING_LAYER_PRIORITY` to prefer
43
+
OpenMP over TBB in `__init__.py`, eliminating fork warnings when Joblib forks
44
+
after Numba parallel execution. Added `llvm-openmp` to conda dependencies.
45
+
-**Fast-path threshold** (`gcate_glm.py`): Raised the effective design-dimension
46
+
ceiling so larger `r` values and wide batch designs correctly use the crispyx path.
47
+
-**Backend toggle** (`gcate_glm.py`): Re-added `_USE_FAST_BACKEND` module flag
48
+
and `_backend_override()` context manager for reliable statsmodels fallback.
49
+
-**Weighted dispersion** (`nb_glm_fast.py`): Dispersion averaging is now
- Default `kwargs_es['max_iters']` reduced from 500 → 50.
74
-
- Default `tolerance` reduced from `1e-3` → `0.0` and a new scale-
75
-
invariant `rel_tol=2e-4` introduced. In practice the relative
76
-
threshold plus `patience=5` reaches the same fixed point as the
77
-
old absolute threshold within ≤ 50 iterations on every tested
78
-
dataset, but external benchmarks that relied on the 500-iteration
79
-
upper bound should pass `kwargs_es_1=dict(max_iters=500)` and
80
-
`kwargs_es_2=dict(max_iters=500)` to recover the pre-v0.0.6 budget.
65
+
- Default `tolerance` reduced from `1e-3` → `0.0`; new scale-invariant
66
+
`rel_tol=2e-4` introduced. To reproduce pre-v0.0.6 behavior, pass
67
+
`kwargs_es_1=dict(max_iters=500)` and `kwargs_es_2=dict(max_iters=500)`.
81
68
82
69
-**⚠️ BREAKING — `LFC()` variance and default `usevar`** (`DR_learner.py`):
83
-
- Default `usevar` changed from `'pooled'` to `'unequal'` (Welch). Callers
84
-
not specifying `usevar=` will see Welch t-statistics instead of pooled
85
-
ones; revert with `LFC(..., usevar='pooled')` if reproducing pre-v0.0.6
86
-
results is required.
87
-
-**`'unequal'` formula corrected**: the variance is now
88
-
`var = s₀²/n₀ + s₁²/n₁` (standard Welch); previously it was
89
-
`(s₀²/n₀ + s₁²/n₁)/2`, a "half-Welch" that under-estimated the
90
-
standard error by √2. Test statistics are therefore ~√2× smaller
91
-
than pre-v0.0.6 for explicit `usevar='unequal'` callers.
92
-
- p-values now use the t-distribution with Welch-Satterthwaite degrees
93
-
of freedom (per gene); the prior version used a Normal approximation.
94
-
95
-
-**`gcate_opt.py`**: `alter_min()` initialisation now uses the auto-dispatch GLM path.
96
-
-**`gcate.py`**: `_check_input()` and `estimate_r()` now use the auto-dispatch dispersion and GLM paths.
97
-
-**`DR_estimation.py`**: `cross_fitting()` now uses the auto-dispatch GLM path, enabling the per-perturbation fast path for multi-treatment LFC estimation.
98
-
-**`DR_learner.py`**: `LFC()` now accepts `backend: str = "auto"` as an explicit parameter; `"fast"` forces crispyx, `"original"` forces statsmodels.
99
-
-**`utils.py`**: `comp_size_factor()` vectorized with `np.nanmean`/`np.nanmedian`.
70
+
- Default `usevar` changed from `'pooled'` to `'unequal'` (Welch). Revert
71
+
with `LFC(..., usevar='pooled')` if reproducing pre-v0.0.6 results.
72
+
-`'unequal'` formula corrected: variance is now `s₀²/n₀ + s₁²/n₁`
73
+
(standard Welch); the prior version used `(s₀²/n₀ + s₁²/n₁)/2`
74
+
("half-Welch"), under-estimating the standard error by √2.
75
+
- p-values now use the t-distribution with Welch-Satterthwaite degrees of
76
+
freedom per gene; the prior version used a Normal approximation.
77
+
78
+
-`alter_min()` initialisation, `_check_input()`, `estimate_r()`, and
79
+
`cross_fitting()` now use the auto-dispatch GLM/dispersion paths.
80
+
-`LFC()` now accepts `backend: str = "auto"` (`"fast"` forces crispyx,
81
+
`"original"` forces statsmodels).
82
+
-`comp_size_factor()` vectorized with `np.nanmean`/`np.nanmedian`.
100
83
101
84
### Performance
102
85
@@ -109,31 +92,19 @@ Benchmarked on Perturb-seq data (n = 2,926 cells, p = 3,221 genes, 29 perturbati
109
92
|**Total**|**419.3 s**|**364.2 s**|**1.2×**|
110
93
111
94
On synthetic data (n = 500, p = 200): 61.5× GLM fit speedup, 7.1× imputation speedup.
95
+
Latent factor recovery: mean canonical correlation 0.998. LFC correlation: 0.856.
112
96
113
-
Latent factor recovery: mean canonical correlation 0.998. LFC correlation: 0.856 (expected difference due to per-perturbation vs. joint model specification).
Three targeted changes reduce end-to-end `gcate_lfc_batch` wall time by **1.48×** on the Replogle tutorial dataset (79,865 cells × 8,563 genes, 200 perturbations, 14 batches):
| Stage 1 ≤3,000-cell mixed subsample | −55 min | tau r=0.992, Jaccard=0.80 |
104
+
| Stage 2 joint fit | −5 min | tau r=0.9994, Jaccard=0.975 |
105
+
|**Combined**|**−70.6 min / 1.48×**|**tau r=0.9994, Jaccard=0.975**|
135
106
136
-
The Stage 1 subsample uses a random draw of up to 3,000 cells from both control and perturbation cells; column scaling and covariate offsets still use all cells. The Stage 2 joint fit is mathematically equivalent to the original per-perturbation loop for one-hot treatment designs (block-diagonal IRLS normal equations).
107
+
Full-run: 217.5 min → 146.9 min (**1.48×** faster); sig pairs −0.2%, perts with ≥1 hit −0.6%.
0 commit comments