Skip to content

Hybrid native backend, protein-aware diffusion, and smarter candidate triage - #14

Closed
cosmic-hydra with Codex wants to merge 2 commits into
mainfrom
codex/upgrade-zane-hybrid-system
Closed

Hybrid native backend, protein-aware diffusion, and smarter candidate triage#14
cosmic-hydra with Codex wants to merge 2 commits into
mainfrom
codex/upgrade-zane-hybrid-system

Conversation

@Codex

@Codex Codex AI commented Apr 19, 2026

Copy link
Copy Markdown
Contributor

Upgrades ZANE into the requested hybrid system: native physics backend, GPU-aware diffusion with protein conditioning, joint interaction modeling, and smarter candidate selection with retrosynthesis-aware filtering.

  • Native physics backend

    • Added C++/CUDA torch extension (compute_energy, compute_forces, run_fep) with CPU fallback and packaging support.
    • Integrated into MD simulator for force/energy/FEP use; exported via physics API.
  • Protein-aware generative modeling

    • Diffusion model now accepts protein/pocket embeddings through cross-attention/FiLM and supports classifier-free guidance.
    • Added protein–ligand interaction head producing affinity and contact maps.
  • Candidate triage and evaluation

    • Selection blends MC-dropout uncertainty, EHVI, Tanimoto diversity, and retrosynthesis feasibility; filters infeasible routes.
    • Added time-based split helper and enrichment factor metric for evaluation.
  • Training stability and regularization

    • Target normalization, label-noise checks, calibration logging, and optional energy regularization hook tied to the native backend.

Example usage (physics-backed energy in training):

trainer = SelfLearningTrainer(
    model,
    energy_regularization_weight=0.1,
    energy_function=lambda batch: compute_energy(batch.pos, reduce=True),
)

Co-authored-by: cosmic-hydra <140935487+cosmic-hydra@users.noreply.github.com>
@cosmic-hydra
cosmic-hydra marked this pull request as ready for review April 19, 2026 10:17
Copilot AI review requested due to automatic review settings April 19, 2026 10:17
@qodo-code-review

Copy link
Copy Markdown

Review Summary by Qodo

Hybrid native backend, protein-aware diffusion, and intelligent candidate triage

✨ Enhancement 🧪 Tests

Grey Divider

Walkthroughs

Description
• Native C++/CUDA physics backend with Lennard-Jones energy, forces, and FEP computation
• Protein-aware diffusion model with cross-attention and FiLM conditioning for ligand generation
• Joint protein–ligand interaction model producing binding affinity and contact maps
• Smarter candidate selection combining MC-dropout uncertainty, EHVI, diversity, and retrosynthesis
  feasibility
• Time-based dataset splitting and MC-dropout uncertainty quantification for predictions
• Training stability improvements: target normalization, calibration tracking, and energy
  regularization
Diagram
flowchart LR
  A["Native Physics Backend<br/>C++/CUDA"] -->|compute_energy<br/>compute_forces<br/>run_fep| B["MD Simulator<br/>& Energy Calc"]
  C["Protein Context<br/>Embeddings"] -->|cross-attention<br/>FiLM| D["Diffusion Model<br/>Ligand Generation"]
  E["Ligand & Residue<br/>Features"] -->|joint encoder| F["Interaction Head<br/>Affinity + Contacts"]
  G["Candidates<br/>with Uncertainty"] -->|EHVI + Diversity<br/>+ Retrosynthesis| H["Ranked Selection<br/>Top-K + Explorers"]
  B -->|coordinates| F
  D -->|generated molecules| G
  F -->|binding scores| H
Loading

Grey Divider

File Changes

1. drug_discovery/native/__init__.py ✨ Enhancement +13/-0

Native backend module initialization and exports

drug_discovery/native/init.py


2. drug_discovery/native/backend.py ✨ Enhancement +193/-0

Torch extension loader with CPU/CUDA fallback implementation

drug_discovery/native/backend.py


3. drug_discovery/native/forcefield.cpp ✨ Enhancement +80/-0

C++/CUDA Lennard-Jones energy, forces, and FEP kernels

drug_discovery/native/forcefield.cpp


View more (13)
4. drug_discovery/models/protein_ligand.py ✨ Enhancement +81/-0

Joint protein–ligand interaction model with affinity and contact heads

drug_discovery/models/protein_ligand.py


5. drug_discovery/models/diffusion_generator.py ✨ Enhancement +128/-3

Protein conditioning via cross-attention and FiLM; classifier-free guidance

drug_discovery/models/diffusion_generator.py


6. drug_discovery/models/__init__.py ✨ Enhancement +3/-0

Export new protein–ligand interaction model classes

drug_discovery/models/init.py


7. drug_discovery/optimization/selection.py ✨ Enhancement +105/-0

Candidate selector blending EHVI, uncertainty, diversity, and feasibility

drug_discovery/optimization/selection.py


8. drug_discovery/optimization/__init__.py ✨ Enhancement +3/-0

Export candidate selection configuration and selector

drug_discovery/optimization/init.py


9. drug_discovery/data/dataset.py ✨ Enhancement +27/-0

Time-based chronological dataset split helper

drug_discovery/data/dataset.py


10. drug_discovery/evaluation/predictor.py ✨ Enhancement +38/-0

MC-dropout uncertainty quantification and enrichment factor metric

drug_discovery/evaluation/predictor.py


11. drug_discovery/physics/__init__.py ✨ Enhancement +4/-0

Export native backend functions from physics module

drug_discovery/physics/init.py


12. drug_discovery/physics/md_simulator.py ✨ Enhancement +120/-24

Integrate native backend for energy/forces; add FEP and coordinate extraction

drug_discovery/physics/md_simulator.py


13. drug_discovery/training/trainer.py ✨ Enhancement +69/-5

Target normalization, calibration tracking, and energy regularization hook

drug_discovery/training/trainer.py


14. drug_discovery/pipeline.py ✨ Enhancement +64/-22

Time-split support, candidate generation with retrosynthesis filtering, energy regularization

drug_discovery/pipeline.py


15. setup.py ⚙️ Configuration changes +2/-0

Include C++ extension source files in package data

setup.py


16. README.md 📝 Documentation +4/-0

Document native backend, protein-aware diffusion, and candidate triage features

README.md


Grey Divider

Qodo Logo

@qodo-code-review

qodo-code-review Bot commented Apr 19, 2026

Copy link
Copy Markdown

Code Review by Qodo

🐞 Bugs (8) 📘 Rule violations (0) 📎 Requirement gaps (0)

Grey Divider


Action required

1. calculate_energy return shape 🐞 Bug ≡ Correctness
Description
EnergyCalculator.calculate_energy() returns plain None in the exception path even when
return_coords=True, so callers that unpack (energy, coords) will crash.
Code

drug_discovery/physics/md_simulator.py[R352-355]

        except Exception as e:
            logger.error(f"Energy calculation error: {e}")
            return None

-    def optimize_geometry(self, smiles: str, max_iters: int = 200) -> tuple[str | None, float | None]:
Evidence
simulate_ligand unconditionally unpacks calculate_energy(..., return_coords=True) into
(initial_energy, coord_tensor). The exception handler in calculate_energy returns None regardless of
return_coords, which causes a TypeError during unpacking when any exception occurs (RDKit embedding,
native backend, etc.).

drug_discovery/physics/md_simulator.py[134-137]
drug_discovery/physics/md_simulator.py[306-355]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

### Issue description
`EnergyCalculator.calculate_energy(..., return_coords=True)` is unpacked by callers, but the exception handler returns a bare `None`, causing an unpacking crash.

### Issue Context
The method already returns `(None, None)` for several early failure cases when `return_coords=True`; the remaining exception handler should follow the same contract.

### Fix Focus Areas
- drug_discovery/physics/md_simulator.py[306-355]

### Suggested fix
- In the `except Exception as e:` block, return `(None, None)` when `return_coords` is `True`, else return `None`.
- Keep logging behavior unchanged.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


2. optimize_geometry return shape 🐞 Bug ≡ Correctness
Description
EnergyCalculator.optimize_geometry() returns a 2-tuple (None, None) on exception even when
return_coords=True, so callers expecting 3 values will crash.
Code

drug_discovery/physics/md_simulator.py[R408-409]

        except Exception as e:
            logger.error(f"Geometry optimization error: {e}")
Evidence
simulate_ligand unpacks optimize_geometry(..., return_coords=True) into three values.
optimize_geometry’s exception handler always returns (None, None), violating the advertised return
type and causing a ValueError/TypeError during unpacking when exceptions occur.

drug_discovery/physics/md_simulator.py[134-137]
drug_discovery/physics/md_simulator.py[356-410]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

### Issue description
`EnergyCalculator.optimize_geometry(..., return_coords=True)` is unpacked as 3 values, but the exception path returns only 2 values.

### Issue Context
The function already returns `(None, None, None)` for some early failures when `return_coords=True`; the exception handler should do the same.

### Fix Focus Areas
- drug_discovery/physics/md_simulator.py[356-410]

### Suggested fix
- In the `except Exception as e:` block, return `(None, None, None)` when `return_coords` is `True`, else return `(None, None)`.
- Ensure type annotation remains accurate.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


3. Energy penalty mixes molecules 🐞 Bug ≡ Correctness
Description
Pipeline energy regularization computes Lennard-Jones energy over concatenated PyG batch.pos,
introducing cross-molecule interactions and an incorrect penalty signal during graph training.
Code

drug_discovery/pipeline.py[R275-282]

+        def _energy_penalty(batch):
+            coords = getattr(batch, "pos", None)
+            if coords is None:
+                return None
+            try:
+                return torch.mean(compute_energy(coords.to(self.device), reduce=False))
+            except Exception:
+                return None
Evidence
For graph training, SelfLearningTrainer passes the PyG Batch object into energy_function. The
pipeline’s _energy_penalty reads batch.pos and calls compute_energy on it directly. compute_energy
treats a 2D (N,3) tensor as one system (adds a batch dim) and computes all-pairs distances via
cdist, so energy includes interactions between atoms belonging to different graphs in the batch.

drug_discovery/pipeline.py[274-293]
drug_discovery/training/trainer.py[105-141]
drug_discovery/native/backend.py[55-71]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

### Issue description
Energy regularization currently treats a PyG minibatch as a single molecule/system, so the LJ penalty is mathematically wrong and mixes unrelated graphs.

### Issue Context
`batch.pos` in PyG is concatenated across graphs; `batch.batch` (or `batch.ptr`) provides the graph membership needed to split coordinates.

### Fix Focus Areas
- drug_discovery/pipeline.py[274-283]
- drug_discovery/native/backend.py[55-71]

### Suggested fix
- In `_energy_penalty`, if `batch` has `batch`/`ptr`, split `coords` into per-graph tensors and compute energy per graph, then average.
 - Example approach: iterate graphs using `batch.ptr` slices; pad to `(B, Nmax, 3)` and call `compute_energy(..., reduce=False)`; or call `compute_energy` per graph in a loop if batches are small.
- If graph membership is unavailable, return `None` (don’t apply penalty).

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


View more (3)
4. MC dropout mutates BatchNorm 🐞 Bug ≡ Correctness
Description
predict_with_uncertainty() sets the whole model to train() for MC dropout, which updates BatchNorm
running stats and can silently change subsequent predictions and training behavior.
Code

drug_discovery/evaluation/predictor.py[R43-54]

+    def predict_with_uncertainty(self, features: torch.Tensor, samples: int = 8) -> tuple[np.ndarray, np.ndarray]:
+        """Predict with MC dropout-based uncertainty."""
+        preds = []
+        was_training = self.model.training
+        self.model.train()
+        with torch.no_grad():
+            for _ in range(max(1, samples)):
+                preds.append(self.model(features.to(self.device)).cpu())
+        self.model.eval()
+        if was_training:
+            self.model.train()
+        stacked = torch.stack(preds, dim=0)
Evidence
PropertyPredictor.predict_with_uncertainty toggles model.train() and runs forward passes. The
primary GNN model uses nn.BatchNorm1d layers in forward; BatchNorm updates running_mean/var in train
mode even under torch.no_grad(), so repeated uncertainty calls drift model state.

drug_discovery/evaluation/predictor.py[43-57]
drug_discovery/models/gnn.py[64-99]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

### Issue description
MC-dropout inference should not update BatchNorm running statistics; current implementation does.

### Issue Context
The GNN uses `nn.BatchNorm1d`, which updates buffers in `train()` mode.

### Fix Focus Areas
- drug_discovery/evaluation/predictor.py[43-57]

### Suggested fix
- Keep the model in `eval()` mode, then enable dropout layers only:
 - Save current mode.
 - Call `self.model.eval()`.
 - Recursively set `Dropout` modules to `train()` (or use a helper that toggles only dropout).
 - Run multiple forward passes.
 - Restore original module modes.
- Alternatively, if you must call `train()`, set all BatchNorm modules back to `eval()` before sampling.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


5. CFG guidance randomly disabled 🐞 Bug ≡ Correctness
Description
Diffusion sampling computes the “unconditional” pass by sometimes reusing the conditioned context,
which frequently makes uncond==cond and reduces classifier-free guidance to a no-op.
Code

drug_discovery/models/diffusion_generator.py[R241-260]

+            if cond_ctx is not None:
+                with torch.no_grad():
+                    if torch.rand(1, device=self.device).item() < self._uncond_prob:
+                        eff_cond_ctx = None
+                        eff_cond_mask = None
+                    else:
+                        eff_cond_ctx = cond_ctx
+                        eff_cond_mask = cond_mask
+                uncond_pos, uncond_atom = self.model(
+                    atom_types,
+                    pos,
+                    edge_index,
+                    t,
+                    batch,
+                    protein_context=eff_cond_ctx,
+                    protein_mask=eff_cond_mask,
+                )
+                eps_pos = uncond_pos + g_scale * (cond_pos - uncond_pos)
+                eps_atom = uncond_atom + g_scale * (cond_atom - uncond_atom)
+            else:
Evidence
In sample(), when cond_ctx is present, the code sometimes sets eff_cond_ctx=cond_ctx for the uncond
forward, so uncond_pos/atom are produced with the same conditioning as cond_pos/atom. In those
steps, eps reduces to the conditional output and guidance is effectively not applied.

drug_discovery/models/diffusion_generator.py[226-262]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

### Issue description
Classifier-free guidance requires both a conditioned and an unconditional prediction at each step; current sampling sometimes computes the 'unconditional' prediction with the conditioning still applied.

### Issue Context
`uncond_dropout_prob` is typically used during training to randomly drop conditioning. During sampling, CFG usually computes:
- `eps_cond = model(x, cond)`
- `eps_uncond = model(x, None)`
- `eps = eps_uncond + s*(eps_cond-eps_uncond)`

### Fix Focus Areas
- drug_discovery/models/diffusion_generator.py[226-262]

### Suggested fix
- In `sample()`, when `cond_ctx` is not None, always compute `uncond_*` with `protein_context=None` (and mask None/zeros) rather than probabilistically.
- Keep `uncond_dropout_prob` for training-time conditioning dropout (if/when training code is added).

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


6. protein_context expand crash 🐞 Bug ☼ Reliability
Description
_prepare_condition() uses Tensor.expand to match num_molecules, which throws when protein_context
has a batch dimension that is neither 1 nor num_molecules.
Code

drug_discovery/models/diffusion_generator.py[R207-211]

+        if ctx.dim() == 2:
+            ctx = ctx.unsqueeze(0).expand(num_molecules, -1, -1)
+        elif ctx.dim() == 3 and ctx.size(0) != num_molecules:
+            ctx = ctx.expand(num_molecules, -1, -1)
+        protein_mask = torch.zeros(ctx.shape[:2], dtype=torch.bool, device=self.device)
Evidence
If a caller provides protein_context shaped (B_ctx, L, D) where B_ctx != 1 and B_ctx !=
num_molecules, .expand(num_molecules, ...) will raise because expand cannot change a non-singleton
dimension.

drug_discovery/models/diffusion_generator.py[203-212]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

### Issue description
`Tensor.expand` only works for singleton dimensions; mismatched batch sizes will crash sampling.

### Issue Context
Users may call `sample(num_molecules=N, protein_context=ctx)` with ctx already batched.

### Fix Focus Areas
- drug_discovery/models/diffusion_generator.py[203-212]

### Suggested fix
- If `ctx.dim()==3`:
 - If `ctx.size(0)==1`, `expand` is fine.
 - If `ctx.size(0)==num_molecules`, keep as-is.
 - Otherwise: raise a clear ValueError explaining the mismatch, or use `repeat` only if semantics are intended.
- Add a unit test or assertion to cover mismatch behavior.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools



Remediation recommended

7. Native FEP None schedule 🐞 Bug ➹ Performance
Description
Python wrapper passes lambda_schedule=None into the C++ extension call, which is not compatible with
the current pybind signature and forces a fallback to the Torch implementation.
Code

drug_discovery/native/backend.py[R171-178]

+    backend = _load_ext()
+    if backend is not None:
+        try:
+            delta_f = backend.run_fep(ligand_coords, protein_coords, lambda_schedule, sigma, epsilon)
+        except Exception as exc:  # pragma: no cover
+            logger.warning("Native FEP failed, falling back to torch: %s", exc)
+            delta_f = _torch_fep(ligand_coords, protein_coords, lambda_schedule, sigma=sigma, epsilon=epsilon)
+    else:
Evidence
The wrapper forwards lambda_schedule directly to backend.run_fep. The C++ binding declares the
argument as torch::Tensor, while the Python API exposes it as Optional; passing None will fail at
the Python↔C++ boundary, triggering the warning + fallback path every time the schedule is omitted.

drug_discovery/native/backend.py[160-180]
drug_discovery/native/forcefield.cpp[47-80]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

### Issue description
Native `run_fep` is effectively unusable with default arguments because `lambda_schedule=None` cannot be passed to a pybind function expecting `torch::Tensor`.

### Issue Context
The C++ implementation already checks for `!lambda_schedule.defined()` / empty schedules, but that logic is unreachable if the binding rejects `None`.

### Fix Focus Areas
- drug_discovery/native/backend.py[160-180]
- drug_discovery/native/forcefield.cpp[47-80]

### Suggested fix
Choose one:
1. **C++/pybind fix (preferred):** change the binding signature to accept an optional tensor (e.g., `c10::optional<torch::Tensor>`), and handle undefined inside C++.
2. **Python wrapper fix:** if `lambda_schedule is None`, pass an empty tensor on the correct device/dtype (so C++ sees `.defined()==true` but `.numel()==0`) and keep the current C++ fallback logic.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


8. EHVI weight/direction bug 🐞 Bug ≡ Correctness
Description
CandidateSelector’s EHVI uses ref - values for all objectives and then adds EHVI with a hardcoded
0.1 factor, which can reward worse (lower) QED and ignores the configured ehvi_weight.
Code

drug_discovery/optimization/selection.py[R39-43]

+def expected_hypervolume_improvement(values: np.ndarray, reference_point: Sequence[float]) -> np.ndarray:
+    ref = np.asarray(reference_point, dtype=np.float32)
+    improvements = np.maximum(ref - values, 0.0)
+    return np.prod(improvements, axis=1)
+
Evidence
The selector’s metrics include qed_score (a maximize objective elsewhere in the codebase). EHVI is
computed as max(ref - values, 0), which increases as qed decreases below ref=1.0, and the
combination step uses +0.1*ehvi rather than CandidateSelectionConfig.ehvi_weight.

drug_discovery/optimization/selection.py[39-43]
drug_discovery/optimization/selection.py[45-96]
drug_discovery/optimization/multi_objective.py[38-44]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

### Issue description
EHVI is computed/weighted inconsistently with maximize-type objectives (e.g., QED) and the config weight is not applied.

### Issue Context
`MultiObjectiveOptimizer` treats `qed_score` as maximize. EHVI should be defined with consistent objective directions (either transform to a minimization space first, or compute improvements per objective accordingly).

### Fix Focus Areas
- drug_discovery/optimization/selection.py[39-43]
- drug_discovery/optimization/selection.py[84-96]

### Suggested fix
- Decide objective directions (min vs max) for each dimension in `metric_values`.
 - For maximize objectives, use `max(values - ref, 0)` (or transform values to minimization form before EHVI).
- Replace `0.1 * ehvi[idx]` with `self.config.ehvi_weight * ehvi[idx]` (and avoid also multiplying `base_metric` by `ehvi_weight` unless that’s intended).

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


Grey Divider

ⓘ The new review experience is currently in Beta. Learn more

Grey Divider

Qodo Logo

Comment on lines 352 to 355
except Exception as e:
logger.error(f"Energy calculation error: {e}")
return None

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

1. Calculate_energy return shape 🐞 Bug ≡ Correctness

EnergyCalculator.calculate_energy() returns plain None in the exception path even when
return_coords=True, so callers that unpack (energy, coords) will crash.
Agent Prompt
### Issue description
`EnergyCalculator.calculate_energy(..., return_coords=True)` is unpacked by callers, but the exception handler returns a bare `None`, causing an unpacking crash.

### Issue Context
The method already returns `(None, None)` for several early failure cases when `return_coords=True`; the remaining exception handler should follow the same contract.

### Fix Focus Areas
- drug_discovery/physics/md_simulator.py[306-355]

### Suggested fix
- In the `except Exception as e:` block, return `(None, None)` when `return_coords` is `True`, else return `None`.
- Keep logging behavior unchanged.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

Comment on lines 408 to 409
except Exception as e:
logger.error(f"Geometry optimization error: {e}")

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

2. Optimize_geometry return shape 🐞 Bug ≡ Correctness

EnergyCalculator.optimize_geometry() returns a 2-tuple (None, None) on exception even when
return_coords=True, so callers expecting 3 values will crash.
Agent Prompt
### Issue description
`EnergyCalculator.optimize_geometry(..., return_coords=True)` is unpacked as 3 values, but the exception path returns only 2 values.

### Issue Context
The function already returns `(None, None, None)` for some early failures when `return_coords=True`; the exception handler should do the same.

### Fix Focus Areas
- drug_discovery/physics/md_simulator.py[356-410]

### Suggested fix
- In the `except Exception as e:` block, return `(None, None, None)` when `return_coords` is `True`, else return `(None, None)`.
- Ensure type annotation remains accurate.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

Comment thread drug_discovery/pipeline.py Outdated
Comment on lines +275 to +282
def _energy_penalty(batch):
coords = getattr(batch, "pos", None)
if coords is None:
return None
try:
return torch.mean(compute_energy(coords.to(self.device), reduce=False))
except Exception:
return None

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

3. Energy penalty mixes molecules 🐞 Bug ≡ Correctness

Pipeline energy regularization computes Lennard-Jones energy over concatenated PyG batch.pos,
introducing cross-molecule interactions and an incorrect penalty signal during graph training.
Agent Prompt
### Issue description
Energy regularization currently treats a PyG minibatch as a single molecule/system, so the LJ penalty is mathematically wrong and mixes unrelated graphs.

### Issue Context
`batch.pos` in PyG is concatenated across graphs; `batch.batch` (or `batch.ptr`) provides the graph membership needed to split coordinates.

### Fix Focus Areas
- drug_discovery/pipeline.py[274-283]
- drug_discovery/native/backend.py[55-71]

### Suggested fix
- In `_energy_penalty`, if `batch` has `batch`/`ptr`, split `coords` into per-graph tensors and compute energy per graph, then average.
  - Example approach: iterate graphs using `batch.ptr` slices; pad to `(B, Nmax, 3)` and call `compute_energy(..., reduce=False)`; or call `compute_energy` per graph in a loop if batches are small.
- If graph membership is unavailable, return `None` (don’t apply penalty).

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

Comment on lines +43 to +54
def predict_with_uncertainty(self, features: torch.Tensor, samples: int = 8) -> tuple[np.ndarray, np.ndarray]:
"""Predict with MC dropout-based uncertainty."""
preds = []
was_training = self.model.training
self.model.train()
with torch.no_grad():
for _ in range(max(1, samples)):
preds.append(self.model(features.to(self.device)).cpu())
self.model.eval()
if was_training:
self.model.train()
stacked = torch.stack(preds, dim=0)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

4. Mc dropout mutates batchnorm 🐞 Bug ≡ Correctness

predict_with_uncertainty() sets the whole model to train() for MC dropout, which updates BatchNorm
running stats and can silently change subsequent predictions and training behavior.
Agent Prompt
### Issue description
MC-dropout inference should not update BatchNorm running statistics; current implementation does.

### Issue Context
The GNN uses `nn.BatchNorm1d`, which updates buffers in `train()` mode.

### Fix Focus Areas
- drug_discovery/evaluation/predictor.py[43-57]

### Suggested fix
- Keep the model in `eval()` mode, then enable dropout layers only:
  - Save current mode.
  - Call `self.model.eval()`.
  - Recursively set `Dropout` modules to `train()` (or use a helper that toggles only dropout).
  - Run multiple forward passes.
  - Restore original module modes.
- Alternatively, if you must call `train()`, set all BatchNorm modules back to `eval()` before sampling.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

Comment on lines +241 to +260
if cond_ctx is not None:
with torch.no_grad():
if torch.rand(1, device=self.device).item() < self._uncond_prob:
eff_cond_ctx = None
eff_cond_mask = None
else:
eff_cond_ctx = cond_ctx
eff_cond_mask = cond_mask
uncond_pos, uncond_atom = self.model(
atom_types,
pos,
edge_index,
t,
batch,
protein_context=eff_cond_ctx,
protein_mask=eff_cond_mask,
)
eps_pos = uncond_pos + g_scale * (cond_pos - uncond_pos)
eps_atom = uncond_atom + g_scale * (cond_atom - uncond_atom)
else:

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

5. Cfg guidance randomly disabled 🐞 Bug ≡ Correctness

Diffusion sampling computes the “unconditional” pass by sometimes reusing the conditioned context,
which frequently makes uncond==cond and reduces classifier-free guidance to a no-op.
Agent Prompt
### Issue description
Classifier-free guidance requires both a conditioned and an unconditional prediction at each step; current sampling sometimes computes the 'unconditional' prediction with the conditioning still applied.

### Issue Context
`uncond_dropout_prob` is typically used during training to randomly drop conditioning. During sampling, CFG usually computes:
- `eps_cond = model(x, cond)`
- `eps_uncond = model(x, None)`
- `eps = eps_uncond + s*(eps_cond-eps_uncond)`

### Fix Focus Areas
- drug_discovery/models/diffusion_generator.py[226-262]

### Suggested fix
- In `sample()`, when `cond_ctx` is not None, always compute `uncond_*` with `protein_context=None` (and mask None/zeros) rather than probabilistically.
- Keep `uncond_dropout_prob` for training-time conditioning dropout (if/when training code is added).

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

Comment on lines +207 to +211
if ctx.dim() == 2:
ctx = ctx.unsqueeze(0).expand(num_molecules, -1, -1)
elif ctx.dim() == 3 and ctx.size(0) != num_molecules:
ctx = ctx.expand(num_molecules, -1, -1)
protein_mask = torch.zeros(ctx.shape[:2], dtype=torch.bool, device=self.device)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

6. Protein_context expand crash 🐞 Bug ☼ Reliability

_prepare_condition() uses Tensor.expand to match num_molecules, which throws when protein_context
has a batch dimension that is neither 1 nor num_molecules.
Agent Prompt
### Issue description
`Tensor.expand` only works for singleton dimensions; mismatched batch sizes will crash sampling.

### Issue Context
Users may call `sample(num_molecules=N, protein_context=ctx)` with ctx already batched.

### Fix Focus Areas
- drug_discovery/models/diffusion_generator.py[203-212]

### Suggested fix
- If `ctx.dim()==3`:
  - If `ctx.size(0)==1`, `expand` is fine.
  - If `ctx.size(0)==num_molecules`, keep as-is.
  - Otherwise: raise a clear ValueError explaining the mismatch, or use `repeat` only if semantics are intended.
- Add a unit test or assertion to cover mismatch behavior.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@codex[agent] resolve conflicts and push changes

Copilot AI left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

This PR upgrades ZANE toward a hybrid system by adding a native (Torch-extension) physics backend, extending diffusion generation with protein conditioning + classifier-free guidance, and introducing smarter candidate triage plus new evaluation/splitting utilities.

Changes:

  • Add a lazily-built Torch C++ extension backend (compute_energy, compute_forces, run_fep) and integrate it into MD simulation and training hooks.
  • Extend diffusion sampling/modeling with protein/pocket conditioning and add a protein–ligand interaction model.
  • Add candidate selection utilities (uncertainty + EHVI + diversity + retrosynthesis filtering), plus time-based splits and an enrichment factor metric.

Reviewed changes

Copilot reviewed 16 out of 16 changed files in this pull request and generated 17 comments.

Show a summary per file
File Description
setup.py Packages native C++ source into the distribution.
drug_discovery/training/trainer.py Adds target normalization, calibration logging, and an energy regularization hook.
drug_discovery/pipeline.py Adds time split support, energy-regularized training wiring, and new candidate triage logic.
drug_discovery/physics/md_simulator.py Integrates native energy/forces/FEP and propagates coordinate tensors.
drug_discovery/physics/init.py Re-exports native physics functions via the physics API.
drug_discovery/optimization/selection.py New candidate selector combining uncertainty, EHVI, and diversity.
drug_discovery/optimization/init.py Exposes new selection utilities from the optimization package.
drug_discovery/native/forcefield.cpp Implements the Torch extension entrypoints (energy/forces/FEP).
drug_discovery/native/backend.py Python wrapper that lazy-loads/compiles the extension with Torch fallbacks.
drug_discovery/native/init.py Public native module exports.
drug_discovery/models/protein_ligand.py New protein–ligand interaction model (affinity + contact map).
drug_discovery/models/diffusion_generator.py Adds protein conditioning and classifier-free guidance to diffusion sampling.
drug_discovery/models/init.py Exposes the new protein–ligand model/config.
drug_discovery/evaluation/predictor.py Adds MC-dropout uncertainty APIs and enrichment factor metric.
drug_discovery/data/dataset.py Adds chronological (time-based) dataset split helper.
README.md Documents new native backend, protein-aware diffusion, interaction head, and triage.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Comment on lines +89 to +95
self.config.ehvi_weight * base_metric + self.config.uncertainty_weight * unc + self.config.diversity_weight * div
)
metric_arr = np.asarray(metric_values, dtype=np.float32)
ehvi = expected_hypervolume_improvement(metric_arr, reference_point=self.config.reference_point)
ranked = []
for idx, score in enumerate(scores):
combined = float(score + 0.1 * ehvi[idx])

Copilot AI Apr 19, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The linear score uses base_metric directly (ehvi_weight * base_metric), but EHVI is later computed from both base_metric and qed_score. This double-counts base_metric while qed_score only affects the much smaller 0.1 * ehvi term; consider making the weighting consistent with the intended multi-objective ranking.

Suggested change
self.config.ehvi_weight * base_metric + self.config.uncertainty_weight * unc + self.config.diversity_weight * div
)
metric_arr = np.asarray(metric_values, dtype=np.float32)
ehvi = expected_hypervolume_improvement(metric_arr, reference_point=self.config.reference_point)
ranked = []
for idx, score in enumerate(scores):
combined = float(score + 0.1 * ehvi[idx])
self.config.uncertainty_weight * unc + self.config.diversity_weight * div
)
metric_arr = np.asarray(metric_values, dtype=np.float32)
ehvi = expected_hypervolume_improvement(metric_arr, reference_point=self.config.reference_point)
ranked = []
for idx, score in enumerate(scores):
combined = float(score + self.config.ehvi_weight * ehvi[idx])

Copilot uses AI. Check for mistakes.
Comment on lines +209 to +210
elif ctx.dim() == 3 and ctx.size(0) != num_molecules:
ctx = ctx.expand(num_molecules, -1, -1)

Copilot AI Apr 19, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ctx.expand(num_molecules, -1, -1) will fail unless ctx.size(0) is 1. If protein_context is already batched (B!=1) and num_molecules differs, this will throw at runtime. Validate that ctx.size(0) is 1 or num_molecules, or tile using repeat when appropriate.

Suggested change
elif ctx.dim() == 3 and ctx.size(0) != num_molecules:
ctx = ctx.expand(num_molecules, -1, -1)
elif ctx.dim() == 3:
if ctx.size(0) == num_molecules:
pass
elif ctx.size(0) == 1:
ctx = ctx.expand(num_molecules, -1, -1)
else:
raise ValueError(
f"protein_context batch size must be 1 or match num_molecules "
f"({num_molecules}), got {ctx.size(0)}"
)
else:
raise ValueError(
f"protein_context must have 2 or 3 dimensions, got {ctx.dim()}"
)

Copilot uses AI. Check for mistakes.
Comment on lines +242 to +256
with torch.no_grad():
if torch.rand(1, device=self.device).item() < self._uncond_prob:
eff_cond_ctx = None
eff_cond_mask = None
else:
eff_cond_ctx = cond_ctx
eff_cond_mask = cond_mask
uncond_pos, uncond_atom = self.model(
atom_types,
pos,
edge_index,
t,
batch,
protein_context=eff_cond_ctx,
protein_mask=eff_cond_mask,

Copilot AI Apr 19, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Classifier-free guidance here is stochastic: the "unconditional" pass sometimes uses the full condition (eff_cond_ctx = cond_ctx), making guidance ineffective and adding unrelated randomness. Guidance should deterministically combine a true unconditional pass (no condition) with the conditional pass; if you want to save compute, make that an explicit option instead of random dropout during sampling.

Suggested change
with torch.no_grad():
if torch.rand(1, device=self.device).item() < self._uncond_prob:
eff_cond_ctx = None
eff_cond_mask = None
else:
eff_cond_ctx = cond_ctx
eff_cond_mask = cond_mask
uncond_pos, uncond_atom = self.model(
atom_types,
pos,
edge_index,
t,
batch,
protein_context=eff_cond_ctx,
protein_mask=eff_cond_mask,
uncond_pos, uncond_atom = self.model(
atom_types,
pos,
edge_index,
t,
batch,
protein_context=None,
protein_mask=None,

Copilot uses AI. Check for mistakes.
Comment on lines +97 to +102
cursor = torch.zeros((batch_size,), device=features.device, dtype=torch.long)
for idx, b in enumerate(batch):
pos = cursor[b].item()
dense[b, pos] = features[idx]
padding_mask[b, pos] = False
cursor[b] += 1

Copilot AI Apr 19, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

_dense_from_batch loops over batch and calls .item() each iteration to index cursor. If batch is on GPU this causes per-element CPU/GPU synchronization and will be a major bottleneck. Prefer a vectorized approach (e.g., torch_geometric.utils.to_dense_batch or scatter-based indexing).

Suggested change
cursor = torch.zeros((batch_size,), device=features.device, dtype=torch.long)
for idx, b in enumerate(batch):
pos = cursor[b].item()
dense[b, pos] = features[idx]
padding_mask[b, pos] = False
cursor[b] += 1
if batch.numel() > 0:
perm = torch.argsort(batch, stable=True)
sorted_batch = batch[perm]
sorted_idx = torch.arange(batch.numel(), device=batch.device, dtype=torch.long)
group_start = torch.ones_like(sorted_batch, dtype=torch.bool)
group_start[1:] = sorted_batch[1:] != sorted_batch[:-1]
start_idx = torch.where(group_start, sorted_idx, torch.zeros_like(sorted_idx))
start_idx = torch.cummax(start_idx, dim=0).values
pos_sorted = sorted_idx - start_idx
pos = torch.empty_like(pos_sorted)
pos[perm] = pos_sorted
dense[batch, pos] = features
padding_mask[batch, pos] = False

Copilot uses AI. Check for mistakes.
Comment thread drug_discovery/evaluation/predictor.py Outdated
Comment on lines +43 to +53
def predict_with_uncertainty(self, features: torch.Tensor, samples: int = 8) -> tuple[np.ndarray, np.ndarray]:
"""Predict with MC dropout-based uncertainty."""
preds = []
was_training = self.model.training
self.model.train()
with torch.no_grad():
for _ in range(max(1, samples)):
preds.append(self.model(features.to(self.device)).cpu())
self.model.eval()
if was_training:
self.model.train()

Copilot AI Apr 19, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Switching the entire model to train() for MC dropout can update BatchNorm (and other training-mode state) even under no_grad, changing the model during inference. Prefer enabling dropout layers only while keeping BatchNorm in eval (or otherwise freezing running stats).

Suggested change
def predict_with_uncertainty(self, features: torch.Tensor, samples: int = 8) -> tuple[np.ndarray, np.ndarray]:
"""Predict with MC dropout-based uncertainty."""
preds = []
was_training = self.model.training
self.model.train()
with torch.no_grad():
for _ in range(max(1, samples)):
preds.append(self.model(features.to(self.device)).cpu())
self.model.eval()
if was_training:
self.model.train()
def _enable_mc_dropout(self) -> list[tuple[torch.nn.Module, bool]]:
"""Enable training mode only for dropout layers and return their previous states."""
dropout_types = (
torch.nn.Dropout,
torch.nn.Dropout1d,
torch.nn.Dropout2d,
torch.nn.Dropout3d,
torch.nn.AlphaDropout,
torch.nn.FeatureAlphaDropout,
)
dropout_states = []
for module in self.model.modules():
if isinstance(module, dropout_types):
dropout_states.append((module, module.training))
module.train()
return dropout_states
def predict_with_uncertainty(self, features: torch.Tensor, samples: int = 8) -> tuple[np.ndarray, np.ndarray]:
"""Predict with MC dropout-based uncertainty."""
preds = []
features = features.to(self.device)
was_training = self.model.training
self.model.eval()
dropout_states = self._enable_mc_dropout()
with torch.no_grad():
for _ in range(max(1, samples)):
preds.append(self.model(features).cpu())
for module, module_was_training in dropout_states:
module.train(module_was_training)
self.model.train(was_training)

Copilot uses AI. Check for mistakes.
Comment thread drug_discovery/native/backend.py Outdated
Comment on lines +172 to +174
if backend is not None:
try:
delta_f = backend.run_fep(ligand_coords, protein_coords, lambda_schedule, sigma, epsilon)

Copilot AI Apr 19, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

lambda_schedule can be None here, but the C++ extension binding expects a torch::Tensor. Passing None will raise and force the torch fallback every time. Convert None to an empty tensor (so the extension can apply its default schedule), or update the binding to accept an optional tensor.

Suggested change
if backend is not None:
try:
delta_f = backend.run_fep(ligand_coords, protein_coords, lambda_schedule, sigma, epsilon)
native_lambda_schedule = (
lambda_schedule
if lambda_schedule is not None
else torch.empty(0, device=ligand_coords.device, dtype=ligand_coords.dtype)
)
if backend is not None:
try:
delta_f = backend.run_fep(
ligand_coords, protein_coords, native_lambda_schedule, sigma, epsilon
)

Copilot uses AI. Check for mistakes.
Comment on lines +110 to +114
batch = []
for b_idx in range(dense.size(0)):
n = int(keep[b_idx].sum().item())
batch.extend([b_idx] * n)
batch_tensor = torch.tensor(batch, device=dense.device, dtype=torch.long)

Copilot AI Apr 19, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

_flatten_from_dense builds batch_tensor via a Python list and then torch.tensor(...). This is slow for large batches and can incur host->device copies. Consider constructing batch_tensor with repeat_interleave using per-batch counts to keep it vectorized and device-friendly.

Suggested change
batch = []
for b_idx in range(dense.size(0)):
n = int(keep[b_idx].sum().item())
batch.extend([b_idx] * n)
batch_tensor = torch.tensor(batch, device=dense.device, dtype=torch.long)
counts = keep.sum(dim=1, dtype=torch.long)
batch_tensor = torch.repeat_interleave(
torch.arange(dense.size(0), device=dense.device, dtype=torch.long),
counts,
)

Copilot uses AI. Check for mistakes.
Comment on lines +43 to +45
def predict_with_uncertainty(self, features: torch.Tensor, samples: int = 8) -> tuple[np.ndarray, np.ndarray]:
"""Predict with MC dropout-based uncertainty."""
preds = []

Copilot AI Apr 19, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

New MC-dropout APIs (predict_with_uncertainty / predict_from_smiles_with_uncertainty) are introduced here, but there are no tests exercising them (e.g., output shapes, restoring model mode, stochasticity). Since tests/test_predictor_comprehensive.py already covers other predictor methods, please add coverage for these new uncertainty methods as well.

Copilot uses AI. Check for mistakes.
Comment thread drug_discovery/evaluation/predictor.py Outdated
Comment on lines +304 to +312
"""Compute enrichment factor at given top fraction."""
if len(y_true) == 0:
return 0.0
n_top = max(1, int(len(y_true) * top_fraction))
order = np.argsort(-y_score.reshape(-1))
top_idx = order[:n_top]
hits_top = float(np.sum(y_true.reshape(-1)[top_idx] > 0))
hit_rate = hits_top / n_top
baseline = float(np.mean(y_true.reshape(-1) > 0)) if np.mean(y_true) != 0 else 1e-8

Copilot AI Apr 19, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

enrichment_factor() introduces new evaluation logic but appears untested. Add unit tests for key edge cases (all-zero labels, all-one labels, ties, very small arrays) to ensure the metric is stable and matches the intended definition.

Suggested change
"""Compute enrichment factor at given top fraction."""
if len(y_true) == 0:
return 0.0
n_top = max(1, int(len(y_true) * top_fraction))
order = np.argsort(-y_score.reshape(-1))
top_idx = order[:n_top]
hits_top = float(np.sum(y_true.reshape(-1)[top_idx] > 0))
hit_rate = hits_top / n_top
baseline = float(np.mean(y_true.reshape(-1) > 0)) if np.mean(y_true) != 0 else 1e-8
"""Compute enrichment factor at given top fraction.
The enrichment factor is defined as the hit rate among the highest-scoring
compounds divided by the overall hit rate in the full set.
Edge cases are handled explicitly to keep the metric stable:
- empty inputs return 0.0
- datasets with no positive labels return 0.0
- ties are ordered deterministically via a stable sort
- very small arrays still evaluate at least one top-ranked item
"""
y_true_flat = np.asarray(y_true).reshape(-1)
y_score_flat = np.asarray(y_score).reshape(-1)
if y_true_flat.size == 0 or y_score_flat.size == 0:
return 0.0
if y_true_flat.size != y_score_flat.size:
raise ValueError("y_true and y_score must have the same number of elements")
n_samples = y_true_flat.size
n_top = min(n_samples, max(1, int(n_samples * top_fraction)))
positives = y_true_flat > 0
baseline = float(np.mean(positives))
if baseline == 0.0:
return 0.0
order = np.argsort(-y_score_flat, kind="mergesort")
top_idx = order[:n_top]
hits_top = float(np.sum(positives[top_idx]))
hit_rate = hits_top / n_top

Copilot uses AI. Check for mistakes.
Comment thread drug_discovery/training/trainer.py Outdated
if self._target_mean is None or self._target_std is None:
self._target_mean = targets.mean()
std = targets.std()
self._target_std = std if std > 1e-6 else torch.tensor(1.0, device=targets.device)

Copilot AI Apr 19, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

std is a 0-d Tensor here, so std > 1e-6 produces a Tensor boolean and the ternary will raise "Boolean value of Tensor is ambiguous" at runtime. Use std.item() for the comparison or clamp std (e.g., std = std.clamp_min(1e-6)).

Suggested change
self._target_std = std if std > 1e-6 else torch.tensor(1.0, device=targets.device)
self._target_std = std.clamp_min(1e-6)

Copilot uses AI. Check for mistakes.
@cosmic-hydra

Copy link
Copy Markdown
Owner

@codex[agent] @claude[agent] please resolve conflicts and commit changes

@cosmic-hydra
cosmic-hydra marked this pull request as draft April 19, 2026 10:27
Co-authored-by: cosmic-hydra <140935487+cosmic-hydra@users.noreply.github.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants