jg/pretty minimal gene exp#10
Conversation
…s for better embeddings
…t and gene_exp scripts
…al features, graph filtering improved, encoder model improved, gene exp embeddings example added, mnist example updated with features
… pyright conflict
|
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
ekiefl
left a comment
There was a problem hiding this comment.
Hey @james-golden-arcadia, nice work overall. Glad to see this is all coming together nicely.
I've branched off your branch and I'm working on "packagizing" your changes (https://github.com/Arcadia-Science/glass-box-umap/tree/ek/pretty-minimal-gene-exp-cleanup). When that's done, we'll merge it into this branch, and then merge this branch into main.
Before I can finish my cleanup, I need some clarity on some things. I've made a few comments below, please help me understand what's going on.
| X: NDArray[np.floating], | ||
| n_neighbors: int = 10, | ||
| metric: str = "cosine", | ||
| metric: str = "euclidean", |
There was a problem hiding this comment.
Interesting! Flagging this to make sure the change is purposeful and not just experimentation
There was a problem hiding this comment.
Should have explained it - euclidean is the default in the original parametric umap and scanpy for gene expression.
But after checking scanpy, cosine distance is preferred for learning the umap embedding directly from the hvg features (with no pca). Maybe that should be automatically set?
| class ResidualMLPBlock(nn.Module): | ||
| """ | ||
| Residual MLP block: x -> Linear -> (Norm) -> Act -> Linear -> (Norm) -> +skip | ||
| Keeps everything bias-free and piecewise-linear (if Norm has no affine). | ||
| """ | ||
|
|
||
| def __init__( | ||
| self, | ||
| dim: int, | ||
| hidden_dim: int | None = None, | ||
| activation: str = "relu", # or "leaky_relu" | ||
| negative_slope: float = 0.01, | ||
| use_norm: bool = True, | ||
| ) -> None: | ||
| super().__init__() | ||
| hidden_dim = hidden_dim or dim | ||
|
|
||
| self.fc1 = nn.Linear(dim, hidden_dim, bias=False) | ||
| self.fc2 = nn.Linear(hidden_dim, dim, bias=False) | ||
|
|
||
| # Norm that does NOT introduce an additive learned offset | ||
| self.norm1 = LayerNormDetached(hidden_dim) if use_norm else nn.Identity() | ||
| self.norm2 = LayerNormDetached(dim) if use_norm else nn.Identity() | ||
|
|
||
| if activation == "leaky_relu": | ||
| # self.act = nn.LeakyReLU(negative_slope=negative_slope) | ||
| # a = negative_slope | ||
|
|
||
| self.act = nn.PReLU() # negative_slope=negative_slope) | ||
| nonlin = "leaky_relu" | ||
| else: | ||
| self.act = nn.ReLU() | ||
| nonlin = "relu" | ||
|
|
||
| # He/Kaiming init for ReLU-family activations | ||
| init.kaiming_uniform_(self.fc1.weight, nonlinearity=nonlin) | ||
| init.kaiming_uniform_(self.fc2.weight, nonlinearity=nonlin) | ||
|
|
||
| def forward(self, x: Tensor) -> Tensor: | ||
| h = self.fc1(x) | ||
| h = self.norm1(h) | ||
| h = self.act(h) | ||
| h = self.fc2(h) | ||
| h = self.norm2(h) | ||
| return x + h | ||
|
|
||
|
|
||
| class DefaultEncoder(nn.Module): | ||
| """ | ||
| Bias-free residual MLP encoder for tabular data. | ||
|
|
||
| - Mean-centering/scaling done outside the module (recommended). | ||
| - All Linear layers use bias=False to preserve your local Jacobian property. | ||
| - Residual blocks greatly improve trainability for depth. | ||
| """ | ||
|
|
||
| def __init__( | ||
| self, | ||
| input_dims: tuple[int, ...], | ||
| n_components: int = 2, | ||
| width: int = 128, | ||
| depth: int = 3, # number of residual blocks | ||
| mlp_ratio: float = 2.0, # inner expansion in each block | ||
| activation: str = "leaky_relu", # or "leaky_relu" | ||
| negative_slope: float = 0.01, | ||
| use_norm: bool = True, | ||
| ) -> None: | ||
| super().__init__() | ||
| in_dim = math.prod(input_dims) | ||
| hidden_dim = int(round(width * mlp_ratio)) | ||
|
|
||
| self.flatten = nn.Flatten() | ||
|
|
||
| # Input projection (bias-free) | ||
| self.in_proj = nn.Linear(in_dim, width, bias=False) | ||
|
|
||
| # Init | ||
| if activation == "leaky_relu": | ||
| a = negative_slope | ||
| nonlin = "leaky_relu" | ||
| else: | ||
| a = 0.0 | ||
| nonlin = "relu" | ||
| init.kaiming_uniform_(self.in_proj.weight, a=a, nonlinearity=nonlin) | ||
|
|
||
| # Residual trunk | ||
| self.blocks = nn.Sequential( | ||
| *[ | ||
| ResidualMLPBlock( | ||
| dim=width, | ||
| hidden_dim=hidden_dim, | ||
| activation=activation, | ||
| negative_slope=negative_slope, | ||
| use_norm=use_norm, | ||
| ) | ||
| for _ in range(depth) | ||
| ] | ||
| ) |
There was a problem hiding this comment.
I'm having trouble understanding the intent regarding activations.
- The activation="leaky_relu" path creates a PReLU, not a LeakyReLU — was the switch to PReLU intentional?
- If PReLU is intentional, should negative_slope be removed since PReLU learns its own slope?
- The kaiming init uses negative_slope but PReLU starts at 0.25 — should these match?
- Is the ReLU branch (activation="relu") actually used anywhere, or can we simplify to always use PReLU?
There was a problem hiding this comment.
The general intent was to switch completely to PReLU (saw better embeddings with it).
- Yes, intentional - should have changed it from leaky.
- You're right that this should be fixed - we can pass negative slope (set to 0.25) for the PReLU initial slope.
- You're absolutely right! Then this negative slope of 0.25 should be passed to kaiming init.
- I think we should always use PReLU.
I will clean this up in a commit.
There was a problem hiding this comment.
How does this differ from GlassBoxUMAP.compute_attributions? The code has lots of gene expression language. Does it have to be or is this generalizable computation?
There was a problem hiding this comment.
This has the mapping from PCs to gene features, which hadn't been implemented yet. I think compute_attributions works properly with direct gene features.
I think it could be generalized to any expression like data (rows and cols), and probably combined in compute_attributions. Will make a fix.
There was a problem hiding this comment.
GlassBoxUMAP.compute_attributions works with both raw features and pca features. If PCA preprocessing was used, then it will project back to raw features.
I think the problem is more that if someone comes with already preprocessed PCA features, then those are the raw features as far as glassbox is concerned
| def prelu_to_leaky(self, model: nn.Module) -> nn.Module: | ||
| """Replace all PReLU modules with LeakyReLU using the learned slopes. | ||
|
|
||
| This is needed for Jacobian computation via ``vmap`` + ``jacrev``, which | ||
| requires stateless activations. | ||
|
|
||
| Args: | ||
| model: The model to convert (not modified in-place). | ||
|
|
||
| Returns: | ||
| A deep copy of the model with PReLU replaced by LeakyReLU. | ||
| """ | ||
| model = copy.deepcopy(model) | ||
| for name, module in model.named_modules(): | ||
| if isinstance(module, nn.PReLU): | ||
| slope = ( | ||
| module.weight.detach().item() | ||
| if module.weight.numel() == 1 | ||
| else module.weight.detach().mean().item() | ||
| ) | ||
| parts = name.split(".") | ||
| parent = model | ||
| for p in parts[:-1]: | ||
| parent = getattr(parent, p) | ||
| # parent = getattr(parent, p) if not p.isdigit() else parent[int(p)] | ||
| # if parts[-1].isdigit(): | ||
| # parent[int(parts[-1])] = nn.LeakyReLU(negative_slope=slope) | ||
| # else: | ||
| setattr(parent, parts[-1], nn.LeakyReLU(negative_slope=slope)) | ||
| return model |
There was a problem hiding this comment.
Why overwrite PReLU? Why not just construct LeakyReLU from the get-go?
There was a problem hiding this comment.
This is a wrinkle of the vmap for instant jacobian. Vmap doesn't work for PReLU, only for leaky relu. But PReLU learns better embeddings. So this hideous chimaera may be a necessary evil.
|
Included in #11 |
PR checklist