From 4da4c35d4e3f59caab29e8e85bcf2d08d5e192dd Mon Sep 17 00:00:00 2001 From: isayev Date: Tue, 4 Aug 2026 10:29:04 -0400 Subject: [PATCH 1/2] test: make the malformed-record test possible to pass The assertion was strengthened from "any valid molecule produced a result" to "both did", which is the right assertion -- a batch that aborts right after the corrupt record satisfies the weaker one. But it was derived from reading the source, and the test needs a real NNP, so it could not be run where it was written. CI's slow tier then failed with {'ethanol'} != {'ethanol', 'propanol'}. The code was never the problem. The input was. RDKit's SDMolSupplier, given "this is not a molecule\n$$$$\n" between two valid records, logs "moving to the beginning of the next molecule" and consumes the following record while resynchronizing: the supplier yields [ethanol, None] and propanol is never handed over at all. No implementation of calc_thermo could have produced a propanol result from that file. The corrupt block is now a well-delimited record -- header, counts line promising two atoms, garbage where the coordinates belong, M END, $$$$ -- so the supplier yields exactly [ethanol, None, propanol]. Measured against this repo's RDKit 2025.09.6, along with a check that iter_thermo_records skips the None and keeps both valid records, so the guarded path is confirmed to fire without loading a model. The docstring now records both supplier behaviors the input depends on, since each one silently makes the test vacuous in a different way. --- tests/test_thermo_reference.py | 44 ++++++++++++++++++++++++---------- 1 file changed, 31 insertions(+), 13 deletions(-) diff --git a/tests/test_thermo_reference.py b/tests/test_thermo_reference.py index 89e4d9d2..411155c1 100644 --- a/tests/test_thermo_reference.py +++ b/tests/test_thermo_reference.py @@ -35,15 +35,23 @@ class TestBatchRobustness: def test_malformed_record_does_not_abort_the_batch(self, job_dir): """A None record between two valid ones must be skipped, not crash. - The corrupt block must sit BETWEEN two valid records, not after the - last one: verified against this repo's RDKit (2025.09.6), a forward - ``SDMolSupplier`` iterator (what ``list(...)`` and ``calc_thermo`` - itself use) only ever surfaces an explicit ``None`` for a corrupt - record when a further valid record follows it in the file. A corrupt - trailing record with nothing after it is silently dropped without - producing a ``None`` at all, so the code path under test (an - unguarded ``None.GetConformer()``) would never fire and the test - would pass today for the wrong reason. + Two properties of RDKit's ``SDMolSupplier`` shape this input, both + measured against this repo's version (2025.09.6) rather than assumed: + + 1. The corrupt block must sit BETWEEN two valid records, not after the + last one. A forward iterator (what ``list(...)`` and ``calc_thermo`` + both use) only surfaces an explicit ``None`` for a corrupt record + when a further valid record follows it; a corrupt trailing record is + dropped silently, so the guarded path would never fire. + + 2. The corrupt block must be a **well-delimited record** -- header, + counts line, garbage where coordinates belong, ``M END``, ``$$$$``. + Loose garbage does not work: on ``"this is not a molecule\\n$$$$\\n"`` + the supplier logs "moving to the beginning of the next molecule" and + **consumes the following record while resynchronizing**, yielding + ``[ethanol, None]`` with propanol never handed over at all. The + assertion below would then be impossible to satisfy for a reason that + has nothing to do with the code under test. """ from Auto3D.ASE.thermo import calc_thermo @@ -52,11 +60,21 @@ def test_malformed_record_does_not_abort_the_batch(self, job_dir): _write_mol(good1, smiles="CCO", name="ethanol") _write_mol(good2, smiles="CCCO", name="propanol") - # Sandwich a deliberately corrupt record between two valid ones. - combined = job_dir / "mixed.sdf" - combined.write_text( - good1.read_text() + "this is not a molecule\n$$$$\n" + good2.read_text() + # A corrupt record that ends cleanly at its own $$$$, so the supplier + # yields exactly [ethanol, None, propanol]. Verified: the counts line + # promises two atoms and the atom block does not deliver them. + corrupt = ( + "corrupt\n" + " RDKit 3D\n" + "\n" + " 2 1 0 0 0 0 0 0 0 0999 V2000\n" + "NOT_A_COORD_LINE\n" + " 1 2 1 0\n" + "M END\n" + "$$$$\n" ) + combined = job_dir / "mixed.sdf" + combined.write_text(good1.read_text() + corrupt + good2.read_text()) out = calc_thermo(str(combined), "AIMNET", use_gpu=False) From 45bfaedb996a0d8e0b3f0cd43f3ece7747dcf912 Mon Sep 17 00:00:00 2001 From: isayev Date: Tue, 4 Aug 2026 11:11:15 -0400 Subject: [PATCH 2/2] refactor!: one conformer filter, gated writers, and drop reasons that say why ## One filter Auto3D carried two duplicate filters with the same criterion but different behavior on malformed input: the energy-clustered one raised KeyError on a record with no E_tot, the legacy all-pairs one tolerated it. That divergence was pinned as a strict xfail three waves ago; it is now fixed and the marker is gone. The suite has no xfailed tests. The order was load-bearing: harden the survivor, pin the oracle with recorded values while both filters still existed, then delete the legacy one. Deleting first would have destroyed the oracle its own tests compared against. A record missing E_tot sorts last, via a (True, 0.0) key against (False, energy), and the 0.0 never participates in a comparison. That is deliberate: E_tot is user-supplied, so a genuine positive energy would otherwise let the placeholder sort first and become top_k's reference conformer. Such records are not dropped -- they reach the within-cluster filter, whose energy guard already treats None as "cannot apply, RMSD alone decides". Since no energy gap can prove such a record is not a duplicate, one of them switches clustering off and compares all pairs, which is what the deleted filter did, so its verdicts are preserved. Inverse assertions included, because a dedup guard that declared everything distinct would satisfy every "these must not merge" test while silently disabling deduplication: two energy-less duplicates still collapse, and an energy-less record still merges with an energy-bearing duplicate. ## Gated writers smiles2smi, encode_ids, decode_ids and select_tautomers take a keyword-only overwrite. Permissive where the caller named the file, restrictive where Auto3D invented it: select_tautomers("/data/results.sdf", k=1) truncated /data/results_top_tautomers.sdf with no gate at all. It now resolves and gates the path before reading its input. Both directions are tested. A gate that refuses correctly but also blocks legitimate writes is not an improvement, so the free paths are asserted too, and confirmed by coverage rather than assumed: each gate's line is executed by the real orchestrator finalize path. ## Drop reasons ranking.py logged "No structure converged" whenever nothing survived, including when every conformer was dropped for stereochemistry. Filters now return a frozen FilterResult(kept, dropped) that rejects unknown reason keys, and that literal survives only when unconverged is the sole reason. The final warning previously hard-coded "the rest were dropped by the connectivity, stereochemistry or energy-window filters" -- a hand-maintained disjunction listing reasons it could not know applied. It now reports the actual tally. ## Selector dispatch The k/window if-elif chain becomes a registry checked against config's SELECTOR_FIELDS at import, so adding a selector to config without wiring it fails at import instead of being silently ignored. Reverting the chain first survived the mutation, since two selectors are behaviorally identical either way. A test now wires a genuine third selector through config, registry and method, which is what makes the mutation red. ## Also Five .smi parsers become one. remove_enantiomers and amend_configuration had no input guard at all: a blank line raised a bare IndexError or a tuple-unpack ValueError, and neither skipped comments. Both now raise a named InputValidationError, and blank and comment lines skip. create_enantiomer's len(keys)==1 special case and its post-loop read of a loop variable -- which worked only by a scoping accident -- become one cursor-based pass, verified bit-identical across 0 to 4 stereocenters including ring and leading positions. ASE/thermo.py used atoms.set_calculator/get_calculator, deprecated in ase 3.22.1, in three places. The output-guard re-export in utils/validation.py is gone now that every call site names utils/output_guard directly. Six docstrings claimed changes landed in 4.0.0 or 4.1.0, versions that will not exist. Verified: 1622 passed, 9 skipped, 0 xfailed, ruff clean. --- src/Auto3D/ASE/geometry.py | 3 +- src/Auto3D/ASE/thermo.py | 17 +- src/Auto3D/SPE.py | 3 +- src/Auto3D/batch_opt/optimization_engine.py | 4 +- src/Auto3D/cli/commands/config.py | 2 +- src/Auto3D/cli/commands/properties.py | 7 +- src/Auto3D/config.py | 2 +- src/Auto3D/constants.py | 2 +- src/Auto3D/filtering.py | 358 ++++++++----- src/Auto3D/id_mapping.py | 43 +- src/Auto3D/ranking.py | 241 +++++++-- src/Auto3D/tautomer.py | 58 +- src/Auto3D/utils/convergence.py | 3 +- src/Auto3D/utils/energy.py | 3 +- src/Auto3D/utils/geometry.py | 4 +- src/Auto3D/utils/output_guard.py | 5 +- src/Auto3D/utils/smi_io.py | 20 +- src/Auto3D/utils/stereochemistry.py | 91 ++-- src/Auto3D/utils/validation.py | 45 +- tests/test_chunk_manager.py | 141 +++++ tests/test_cli_errors.py | 2 +- tests/test_cli_property_commands.py | 2 +- tests/test_durability.py | 10 +- tests/test_filter_unique.py | 268 ---------- tests/test_filtering.py | 407 ++++++++++++-- tests/test_id_mapping.py | 10 +- tests/test_output_overwrite_gates.py | 259 +++++++++ tests/test_padding_invariance.py | 11 +- tests/test_ranking.py | 562 +++++++++++++++++--- tests/test_stereo_postopt.py | 24 +- tests/test_utils_stereochemistry.py | 84 +++ tests/test_utils_validation.py | 108 ---- 32 files changed, 1959 insertions(+), 840 deletions(-) delete mode 100644 tests/test_filter_unique.py create mode 100644 tests/test_output_overwrite_gates.py diff --git a/src/Auto3D/ASE/geometry.py b/src/Auto3D/ASE/geometry.py index 8ab14fbc..c09e938a 100644 --- a/src/Auto3D/ASE/geometry.py +++ b/src/Auto3D/ASE/geometry.py @@ -20,11 +20,10 @@ from Auto3D.torch_config import TorchConfig, configure_torch from Auto3D.utils.atomic_io import atomic_write_path from Auto3D.utils.energy import E_TOT_HARTREE_PROP, E_TOT_PROP +from Auto3D.utils.output_guard import check_output_not_input, check_output_overwrite from Auto3D.utils.validation import ( check_engine_supports_molecules, check_gpu_requested, - check_output_not_input, - check_output_overwrite, ) __all__ = ["opt_geometry"] diff --git a/src/Auto3D/ASE/thermo.py b/src/Auto3D/ASE/thermo.py index c535dbe6..9b0cc3dc 100644 --- a/src/Auto3D/ASE/thermo.py +++ b/src/Auto3D/ASE/thermo.py @@ -42,11 +42,10 @@ from Auto3D.torch_config import TorchConfig, configure_torch from Auto3D.utils.energy import hartree2ev from Auto3D.utils.logging_config import get_logger +from Auto3D.utils.output_guard import check_output_not_input, check_output_overwrite from Auto3D.utils.validation import ( check_engine_supports_molecules, check_gpu_requested, - check_output_not_input, - check_output_overwrite, ) __all__ = ["calc_thermo"] @@ -686,7 +685,11 @@ def vib_hessian(mol: Chem.Mol, ase_calculator, model, # mass weighting, and the rotational partition function silently disagree # for isotopically labeled input. atoms = mol2atoms(mol, positions=positions) - atoms.set_calculator(ase_calculator) + # atoms.set_calculator() is deprecated since ase 3.22.1 in favor of the + # `.calc` attribute (Minor 6); `pyproject.toml` pins no ase upper bound + # and globally ignores DeprecationWarning, so removal would otherwise + # land as a silent-until-runtime AttributeError with no advance warning. + atoms.calc = ase_calculator charge = rdmolops.GetFormalCharge(mol) # get the Hessian @@ -1194,7 +1197,9 @@ def do_mol_thermo(mol: Chem.Mol, # them via the explicit `positions=` argument, not from mol's conformer), # so nothing here depends on mol's conformer being in sync yet. coord = atoms.get_positions() - vib = vib_hessian(mol, atoms.get_calculator(), model, device, + # atoms.get_calculator() is deprecated since ase 3.22.1; use `.calc` + # (Minor 6, same rationale as the set_calculator() call above). + vib = vib_hessian(mol, atoms.calc, model, device, model_name=model_name, positions=coord) e = atoms.get_potential_energy() geometry = _detect_geometry(atoms) @@ -1691,7 +1696,9 @@ def calc_thermo(path: str, model_name: str, mol_info_func=None, atoms = mol2atoms(mol) calculator.set_charge(charge) - atoms.set_calculator(calculator) + # atoms.set_calculator() is deprecated since ase 3.22.1; use `.calc` + # (Minor 6, same rationale as vib_hessian's call above). + atoms.calc = calculator if mol_info_func is None: idx = mol.GetProp("_Name").strip() diff --git a/src/Auto3D/SPE.py b/src/Auto3D/SPE.py index 09ad68ca..043f5fe6 100644 --- a/src/Auto3D/SPE.py +++ b/src/Auto3D/SPE.py @@ -13,11 +13,10 @@ from Auto3D.models.preflight import resolve_engine_name from Auto3D.torch_config import TorchConfig, configure_torch from Auto3D.utils.logging_config import get_logger +from Auto3D.utils.output_guard import check_output_not_input, check_output_overwrite from Auto3D.utils.validation import ( check_engine_supports_molecules, check_gpu_requested, - check_output_not_input, - check_output_overwrite, ) logger = get_logger(__name__) diff --git a/src/Auto3D/batch_opt/optimization_engine.py b/src/Auto3D/batch_opt/optimization_engine.py index 289a8f2c..ec6f8426 100644 --- a/src/Auto3D/batch_opt/optimization_engine.py +++ b/src/Auto3D/batch_opt/optimization_engine.py @@ -272,7 +272,7 @@ def _step_active_subset( not_oscillating = oscillating_count < patience # Combine the convergence criteria. An `& ~energy_converged` term stood - # here until 4.0.0; `energy_converged` required `fmax < opttol` while + # here until 3.0.0; `energy_converged` required `fmax < opttol` while # `not_converged_post1` is `fmax > opttol`, so the term was the identity # of `&` wherever it was consulted and false-dominated elsewhere -- it # could never change an outcome, including at the `fmax == opttol` @@ -373,7 +373,7 @@ def n_steps( 2. Oscillation detection: drops structures that don't improve for patience steps There is deliberately no energy-stability criterion. One existed until - 4.0.0 but could never fire: it required ``fmax < opttol``, which is + 3.0.0 but could never fire: it required ``fmax < opttol``, which is exactly the condition under which the force criterion has already stopped the structure, so the term was the identity of ``&`` at every element (see ``test_convergence_outcome_never_depends_on_energy_stability``). diff --git a/src/Auto3D/cli/commands/config.py b/src/Auto3D/cli/commands/config.py index 59dd350c..52414a58 100644 --- a/src/Auto3D/cli/commands/config.py +++ b/src/Auto3D/cli/commands/config.py @@ -119,7 +119,7 @@ def execute_config_init( Both refusals below are configuration problems, so both leave through ``handle_error`` at exit **2**. The overwrite refusal in particular had to move: the CLI's own ``-o`` overwrite gate - (``utils.validation.check_output_overwrite``) raises ``ConfigurationError`` + (``utils.output_guard.check_output_overwrite``) raises ``ConfigurationError`` and exits 2, and the CHANGELOG described the two as printing "the same message" -- while this one hard-coded ``SystemExit(1)``, so a script branching on 2 saw one of them and not the other. diff --git a/src/Auto3D/cli/commands/properties.py b/src/Auto3D/cli/commands/properties.py index 24a73a60..ecac1ab9 100644 --- a/src/Auto3D/cli/commands/properties.py +++ b/src/Auto3D/cli/commands/properties.py @@ -18,11 +18,8 @@ from Auto3D.cli.errors import handle_error from Auto3D.exceptions import ConfigurationError, DependencyError from Auto3D.models.preflight import resolve_engine_name -from Auto3D.utils.validation import ( - check_gpu_requested, - check_output_not_input, - check_output_overwrite, -) +from Auto3D.utils.output_guard import check_output_not_input, check_output_overwrite +from Auto3D.utils.validation import check_gpu_requested # Engine names offered for shell completion. Free-form registry names and custom # model paths are also accepted -- each command below validates them with diff --git a/src/Auto3D/config.py b/src/Auto3D/config.py index 189eefef..22e2f51c 100644 --- a/src/Auto3D/config.py +++ b/src/Auto3D/config.py @@ -484,7 +484,7 @@ class OptimizationConfig: memory but may be faster.""" # There is deliberately no energy_tol/energy_patience here. Both existed - # until 4.0.0 and reached an optimizer criterion that could never fire + # until 3.0.0 and reached an optimizer criterion that could never fire # (audit M1), so they were knobs that changed nothing. Removed rather than # kept as inert configuration. diff --git a/src/Auto3D/constants.py b/src/Auto3D/constants.py index 5047c3b0..698a60f3 100644 --- a/src/Auto3D/constants.py +++ b/src/Auto3D/constants.py @@ -71,7 +71,7 @@ DEFAULT_PATIENCE = 250 # Steps before dropping oscillating conformer DEFAULT_BATCHSIZE_ATOMS = 1024 # Atoms per batch for GPU optimization DEFAULT_ENERGY_CLUSTER_WINDOW = 0.1 # eV, for RMSD clustering -# DEFAULT_ENERGY_TOL / DEFAULT_ENERGY_PATIENCE were removed in 4.0.0 along with +# DEFAULT_ENERGY_TOL / DEFAULT_ENERGY_PATIENCE were removed in 3.0.0 along with # the optimizer's energy-stability criterion, which could never fire (audit M1). # Tuning the tolerance -- for fp32 noise or anything else -- changed nothing. DEFAULT_RANDOM_SEED = 42 # Default random seed for reproducibility diff --git a/src/Auto3D/filtering.py b/src/Auto3D/filtering.py index b4af3b98..7bbe02db 100644 --- a/src/Auto3D/filtering.py +++ b/src/Auto3D/filtering.py @@ -1,7 +1,26 @@ #!/usr/bin/env python -"""Optimized conformer filtering with hierarchical RMSD comparison.""" +"""The one conformer filter: duplicate removal with hierarchical RMSD comparison. + +Auto3D carried two implementations of this until 3.0.0 -- this energy-clustered +one and a legacy all-pairs ``filter_unique``, selected by +``ConformerRanker(use_optimized_filtering=...)``. They applied the identical +duplicate criterion and were kept side by side so each could act as the other's +oracle, but a boolean kwarg that swaps one filter for another is two behaviors +to keep in step, and they had already drifted on malformed input: the legacy one +tolerated a record with no usable ``E_tot`` while this one raised ``KeyError`` +from its sort key. The survivor tolerates it (see :func:`_energy_sort_key`), and +the flag and the duplicate implementation are gone. + +The filter reports **why** it dropped things, not just what survived +(:class:`FilterResult`). ``ranking`` used to log "No structure converged" for a +species whose conformers were all dropped for *stereochemistry* -- a message +that sent the reader to the optimizer settings for a problem in the input's +stereo definitions. +""" from __future__ import annotations +from dataclasses import dataclass + from rdkit import Chem from rdkit.Chem import rdMolAlign @@ -12,67 +31,187 @@ ) from Auto3D.utils.connectivity import check_connectivity from Auto3D.utils.convergence import converged_or_unfiltered -from Auto3D.utils.energy import e_tot_ev, try_e_tot_ev +from Auto3D.utils.energy import try_e_tot_ev from Auto3D.utils.stereo_check import species_key, stereo_preserved +__all__ = [ + "DROP_REASONS", + "FilterResult", + "filter_conformers", + "filter_unique_optimized", +] + +#: Every reason a conformer can leave the selected set, in report order. +#: +#: The authoritative vocabulary: :class:`FilterResult` refuses a count keyed by +#: anything else, so a producer that invents a reason fails at construction +#: rather than contributing a silently unlabeled drop to a diagnostic. +#: ``energy_window`` is produced by ``ranking.ConformerRanker.top_window`` +#: rather than by this module -- it is a selection criterion, not a validity +#: one -- but it belongs in the same vocabulary because it reaches the user +#: through the same message. +DROP_REASONS: tuple[str, ...] = ( + "unparsed", + "unconverged", + "stereochemistry", + "connectivity", + "duplicate", + "energy_window", +) -def filter_unique_optimized( +#: Human-readable phrase per reason, used by :meth:`FilterResult.summary`. +_REASON_PHRASES: dict[str, str] = { + "unparsed": "unparseable by RDKit", + "unconverged": "marked Converged=false", + "stereochemistry": "changed stereochemistry during optimization", + "connectivity": "have broken or newly formed bonds", + "duplicate": "duplicates of a kept conformer", + "energy_window": "outside the energy window", +} + + +@dataclass(frozen=True) +class FilterResult: + """What survived filtering, and a count of what did not, by reason. + + Deliberately tiny: a list and a ``{reason: count}`` dict. The alternative + considered -- attaching a reason to each dropped molecule -- would keep + every rejected conformer alive for the duration of a chunk, which is the + memory the filter exists to release. + + Args: + kept: Surviving molecules, sorted by energy (lowest first); records + with no usable energy sort last. + dropped: Count per reason. Keys must come from :data:`DROP_REASONS`; + reasons that did not fire may be omitted or present as 0. + + Raises: + ValueError: ``dropped`` carries a key outside :data:`DROP_REASONS`. + """ + + kept: list[Chem.Mol] + dropped: dict[str, int] + + def __post_init__(self) -> None: + unknown = sorted(set(self.dropped) - set(DROP_REASONS)) + if unknown: + raise ValueError( + f"unknown filter drop reason(s) {unknown}; expected one of " + f"{list(DROP_REASONS)}" + ) + + @property + def reasons(self) -> tuple[str, ...]: + """The reasons that actually fired, in :data:`DROP_REASONS` order.""" + return tuple(r for r in DROP_REASONS if self.dropped.get(r)) + + def summary(self) -> str: + """One phrase per reason that fired, e.g. ``"2 marked Converged=false"``. + + Empty string when nothing was dropped, so a caller can test it for + truth rather than special-casing a placeholder. + """ + return ", ".join( + f"{self.dropped[reason]} {_REASON_PHRASES[reason]}" + for reason in self.reasons + ) + + +def _energy_sort_key(mol: Chem.Mol) -> tuple[bool, float]: + """Sort key that tolerates a record with no usable ``E_tot``. + + Returns ``(False, energy_ev)`` for a record that has one and + ``(True, 0.0)`` for a record that does not, so the energy-less records land + **after** every record that has an energy, whatever their values. + + The 0.0 is a tie-break placeholder, never an energy: reading a missing + ``E_tot`` as 0.0 on its own would sort a garbage record ahead of every + genuine structure (real ``E_tot`` values are large and negative), making it + the reference conformer ``E_rel`` is measured from and the single structure + a ``k=1`` request returns. Among themselves, energy-less records keep their + input order (``list.sort`` is stable), which is the only ordering there is + any evidence for. + """ + energy = try_e_tot_ev(mol) + if energy is None: + return (True, 0.0) + return (False, energy) + + +def filter_conformers( mols: list[Chem.Mol], + *, rmsd_threshold: float = DEFAULT_RMSD_THRESHOLD, energy_cluster_window: float = DEFAULT_ENERGY_CLUSTER_WINDOW, -) -> list[Chem.Mol]: +) -> FilterResult: """Remove duplicate conformers, skipping RMSD comparisons that cannot match. - Sorts by energy and only RMSD-compares molecules close enough in energy to be - duplicates at all, which avoids the O(n^2) comparisons of the legacy - :func:`filter_unique` below **without changing which molecules - survive** -- the partitioning is chosen so that no duplicate pair can be - separated by it. See the comment on the split rule below for why that - holds; it did not hold before 4.0.0. + Sorts by energy and only RMSD-compares molecules close enough in energy to + be duplicates at all, which avoids comparing all pairs **without changing + which molecules survive** -- the partitioning is chosen so that no duplicate + pair can be separated by it. See the comment on the split rule below for why + that holds; it did not hold before 3.0.0. A pair counts as a duplicate only when all three of these agree: the two are the same compound (:func:`Auto3D.utils.stereo_check.species_key`), their heavy-atom RMSD is under ``rmsd_threshold``, and their energies agree within - ``DEFAULT_DUPLICATE_ENERGY_TOL``. + ``DEFAULT_DUPLICATE_ENERGY_TOL`` (or at least one of them has no usable + energy, in which case that term cannot apply and RMSD alone decides). Args: mols: List of RDKit Mol objects with 'E_tot' (Hartree) and, optionally, 'Converged' properties. A record whose 'Converged' property is explicitly false is dropped; a record without the property is kept (it is not filtered on convergence). Records marked - 'Stereo_changed' are excluded. + 'Stereo_changed' are excluded. ``None`` entries -- what + ``SDMolSupplier`` yields for a record RDKit cannot parse -- are + counted and skipped. rmsd_threshold: RMSD threshold for considering structures similar (Angstrom). energy_cluster_window: Energy width (eV) below which molecules are compared to each other. A performance knob only: values below the duplicate energy tolerance cannot shrink the comparison set, because that tolerance is the floor at which a pair can still be a duplicate. The stored Hartree energies are converted to eV on read - (``Auto3D.utils.energy.e_tot_ev``), so the unit is as documented. + (``Auto3D.utils.energy``), so the unit is as documented. Returns: - List of unique molecules, sorted by energy (lowest first). + A :class:`FilterResult` whose ``kept`` list is sorted by energy (lowest + first) and whose ``dropped`` counts say why the rest are missing. """ + dropped: dict[str, int] = {} + + def _drop(reason: str) -> None: + dropped[reason] = dropped.get(reason, 0) + 1 + # Filter converged structures with valid connectivity. A record with no # 'Converged' property is not filtered on convergence (see # Auto3D.utils.convergence): treating its absence as failure deleted every # record of any SDF batchopt did not write. - valid_mols = [] + # + # Checked in this order, one reason attributed per record, so a structure + # that fails several is reported under the first -- the same short-circuit + # order the single `and` chain here used to have, hence the same verdicts. + valid_mols: list[Chem.Mol] = [] for mol in mols: if mol is None: - continue - if ( - converged_or_unfiltered(mol) - and stereo_preserved(mol) - and check_connectivity(mol) - ): + _drop("unparsed") + elif not converged_or_unfiltered(mol): + _drop("unconverged") + elif not stereo_preserved(mol): + _drop("stereochemistry") + elif not check_connectivity(mol): + _drop("connectivity") + else: valid_mols.append(mol) if not valid_mols: - return [] + return FilterResult(kept=[], dropped=dropped) # Sort by energy. E_tot is stored in Hartree; energy_cluster_window and # the duplicate tolerance below are both in eV, so convert on read. - valid_mols.sort(key=e_tot_ev) + # Records with no usable energy sort last -- see _energy_sort_key. + valid_mols.sort(key=_energy_sort_key) + energies = [try_e_tot_ev(mol) for mol in valid_mols] # Partition the energy-sorted list into runs, and only RMSD-compare within a # run. Where a run may end is a correctness question, not a tuning one. @@ -97,22 +236,34 @@ def filter_unique_optimized( # splitting on LARGER gaps is always safe (it only merges runs and compares # more pairs). It therefore stays a performance knob and can no longer become # a correctness hole. Runs are no longer width-bounded, so a dense energy - # ladder degrades to the O(n^2) of the legacy `filter_unique` -- the price of - # the guarantee, on conformer counts that are tens per species. - split_gap = max(DEFAULT_DUPLICATE_ENERGY_TOL, energy_cluster_window) - clusters: list[list[Chem.Mol]] = [] - current_cluster: list[Chem.Mol] = [valid_mols[0]] - previous_e = e_tot_ev(valid_mols[0]) - - for mol in valid_mols[1:]: - e = e_tot_ev(mol) - if e - previous_e <= split_gap: - current_cluster.append(mol) - else: - clusters.append(current_cluster) - current_cluster = [mol] - previous_e = e - clusters.append(current_cluster) + # ladder degrades to comparing all pairs -- the price of the guarantee, on + # conformer counts that are tens per species. + # + # That entire argument rests on every record HAVING an energy. As soon as one + # does not, no gap proves anything about it: `_filter_within_cluster` falls + # back to RMSD-only for a pair where either side's energy is missing, so such + # a record can be a duplicate of any same-species structure at any energy. + # The honest response is to stop partitioning and compare all pairs, which is + # exactly what the legacy all-pairs filter did with this input -- so the + # survivor keeps its verdicts. Malformed input is rare and reaches here only + # through a direct API call (`ConformerRanker` refuses a record with no + # 'E_tot' up front), so the quadratic cost is paid where it is warranted. + if any(energy is None for energy in energies): + clusters: list[list[Chem.Mol]] = [valid_mols] + else: + split_gap = max(DEFAULT_DUPLICATE_ENERGY_TOL, energy_cluster_window) + clusters = [] + current_cluster: list[Chem.Mol] = [valid_mols[0]] + previous_e = energies[0] + + for mol, e in zip(valid_mols[1:], energies[1:], strict=True): + if e - previous_e <= split_gap: + current_cluster.append(mol) + else: + clusters.append(current_cluster) + current_cluster = [mol] + previous_e = e + clusters.append(current_cluster) # Filter unique within each cluster unique_mols: list[Chem.Mol] = [] @@ -120,7 +271,32 @@ def filter_unique_optimized( unique_in_cluster = _filter_within_cluster(cluster, rmsd_threshold) unique_mols.extend(unique_in_cluster) - return unique_mols + n_duplicates = len(valid_mols) - len(unique_mols) + if n_duplicates: + dropped["duplicate"] = n_duplicates + return FilterResult(kept=unique_mols, dropped=dropped) + + +def filter_unique_optimized( + mols: list[Chem.Mol], + rmsd_threshold: float = DEFAULT_RMSD_THRESHOLD, + energy_cluster_window: float = DEFAULT_ENERGY_CLUSTER_WINDOW, +) -> list[Chem.Mol]: + """The surviving molecules from :func:`filter_conformers`, nothing else. + + Kept as the public name it has always been, for callers that want the list + and not the drop counts. Prefer :func:`filter_conformers` when the reason + something is missing has to reach a user. + + Returns: + List of unique molecules, sorted by energy (lowest first); records with + no usable energy sort last. + """ + return filter_conformers( + mols, + rmsd_threshold=rmsd_threshold, + energy_cluster_window=energy_cluster_window, + ).kept def _filter_within_cluster( @@ -206,105 +382,3 @@ def _mol_energy(mol: Chem.Mol) -> float | None: """ return try_e_tot_ev(mol) - -def filter_unique(mols: list[Chem.Mol], crit: float = DEFAULT_RMSD_THRESHOLD) -> list[Chem.Mol]: - """Remove structures that are very similar and remove unconverged structures. - - The legacy all-pairs filter, kept beside :func:`filter_unique_optimized` - (which supersedes it) because ``ConformerRanker(use_optimized_filtering=False)`` - still selects it and because it is the oracle the optimized filter is - compared against. - - This function filters a list of molecules to keep only unique, converged structures. - It first removes unconverged structures and those with invalid connectivity, - then removes similar structures based on RMSD comparison. - - Args: - mols: List of RDKit molecule objects, optionally carrying a 'Converged' - property. A record whose 'Converged' is explicitly false is - dropped; a record without the property is kept (not filtered on - convergence). Records marked 'Stereo_changed' are excluded. - crit: RMSD threshold for considering two structures as identical. - Structures with RMSD below this value are considered duplicates. - Defaults to DEFAULT_RMSD_THRESHOLD (0.3 Angstroms). - - Returns: - List of unique, converged molecules with valid connectivity. - - Example: - >>> from rdkit import Chem - >>> from rdkit.Chem import AllChem - >>> mol = Chem.MolFromSmiles("CCO") - >>> mol = Chem.AddHs(mol) - >>> AllChem.EmbedMolecule(mol, randomSeed=42) - 0 - >>> mol.SetProp("Converged", "true") - >>> filter_unique([mol], crit=0.3) # Returns list with 1 molecule - [...] - """ - # Remove structures that explicitly failed to converge. A record with no - # 'Converged' property is NOT filtered on convergence -- see - # Auto3D.utils.convergence for why absence is not failure. - mols_: list[Chem.Mol] = [] - for mol in mols: - convergence_flag = converged_or_unfiltered(mol) - has_valid_bonds = check_connectivity(mol) - if convergence_flag and has_valid_bonds and stereo_preserved(mol): - mols_.append(mol) - mols = mols_ - - # Remove similar structures. Strip Hs once per molecule (O(n)) instead of on - # both sides of every comparison (O(n^2)); GetBestRMS on no-H forms is - # symmetric so results are unchanged. The ORIGINAL (H-explicit) mols are - # returned; no-H forms are comparison-only. - # - # Heavy-atom RMSD alone collapses conformers that differ only in an O-H / N-H - # rotor orientation. Guard with an energy check: a pair counts as duplicate - # only when the RMSD is below ``crit`` AND the energies agree within - # DEFAULT_DUPLICATE_ENERGY_TOL (eV; 'E_tot' is stored in Hartree and is - # converted on read by Auto3D.utils.energy). Mols without a usable 'E_tot' - # fall back to RMSD-only (energy guard cannot apply). - unique_mols: list[Chem.Mol] = [] - unique_noH: list[Chem.Mol] = [] - unique_energies: list[float | None] = [] - unique_species: list[str] = [] - for mol_i in mols: - mol_i_noH = Chem.RemoveHs(mol_i) - # E_tot is stored in Hartree; DEFAULT_DUPLICATE_ENERGY_TOL is in eV. - e_i: float | None = try_e_tot_ev(mol_i) - species_i = species_key(mol_i) - unique = True - for mol_j_noH, e_j, species_j in zip( - unique_noH, unique_energies, unique_species, strict=True - ): - # Two different compounds are never duplicates of each other, however - # close their geometries. All stereoisomers of one input share a - # ranking group, and two ring diastereomers can sit below the default - # 0.3 A threshold, so without this the RMSD test could delete one of - # them (Auto3D.utils.stereo_check.species_key). Checked before the - # RMSD call it makes unnecessary. - if species_i != species_j: - continue - try: - # temporary bug fix for https://github.com/rdkit/rdkit/issues/6826 - # removing Hs speeds up the calculation - rmsd = rdMolAlign.GetBestRMS(mol_i_noH, mol_j_noH) - except RuntimeError: - # Incomparable pair: treat as distinct (not a duplicate) so the - # conformer is kept. Using 0 would make it look like a perfect - # duplicate and drop a genuinely distinct structure. - rmsd = float("inf") - energy_close = ( - e_i is None - or e_j is None - or abs(e_i - e_j) < DEFAULT_DUPLICATE_ENERGY_TOL - ) - if rmsd < crit and energy_close: - unique = False - break - if unique: - unique_mols.append(mol_i) - unique_noH.append(mol_i_noH) - unique_energies.append(e_i) - unique_species.append(species_i) - return unique_mols diff --git a/src/Auto3D/id_mapping.py b/src/Auto3D/id_mapping.py index 7b2c45e4..028cf96d 100644 --- a/src/Auto3D/id_mapping.py +++ b/src/Auto3D/id_mapping.py @@ -16,15 +16,19 @@ from rdkit import Chem -from Auto3D.exceptions import ConfigurationError, InputValidationError +from Auto3D.exceptions import InputValidationError from Auto3D.utils.logging_config import get_logger +from Auto3D.utils.output_guard import check_output_overwrite from Auto3D.utils.smi_io import iter_smi_records logger = get_logger(__name__) def encode_ids( - path: str, out_dir: str | os.PathLike[str] | None = None + path: str, + out_dir: str | os.PathLike[str] | None = None, + *, + overwrite: bool = False, ) -> tuple[str, dict[str, int]]: """Encode molecule IDs to numeric indices. @@ -46,6 +50,13 @@ def encode_ids( path: Path to the input .smi or .sdf file. out_dir: Directory to write the encoded file into. Defaults to the input file's own directory. + overwrite: Allow replacing an existing file at the encoded path. + Keyword-only, defaults to False -- Auto3D derives this name, so it + belongs to whoever already owns it. This used to refuse + *unconditionally*, with no way for a caller who meant it to say so; + the keyword and its default now match ``decode_ids`` and + ``tautomer.select_tautomers``, and all three raise through the same + ``utils.output_guard.check_output_overwrite``. Returns: Tuple containing: @@ -54,7 +65,8 @@ def encode_ids( Raises: ValueError: If the input file is neither .smi nor .sdf format. - ConfigurationError: If a file already exists at the encoded path. + ConfigurationError: If a file already exists at the encoded path and + `overwrite` is False. InputValidationError: If a molecule has a missing/blank ID or a duplicate ID is encountered. @@ -72,13 +84,7 @@ def encode_ids( directory = Path(out_dir) if out_dir is not None else path_obj.parent new_path = directory / f"{path_obj.stem}_encoded.{extension}" - if new_path.exists(): - raise ConfigurationError( - f"encode_ids would overwrite the existing file {new_path}. " - "Auto3D writes its encoded copy of the input there; move or " - "rename that file, or pass out_dir to write the encoded copy " - "somewhere else." - ) + check_output_overwrite(new_path, overwrite) if extension == "smi": new_data: list[str] = [] @@ -129,7 +135,9 @@ def encode_ids( return str(new_path), mapping -def decode_ids(path: str, mapping: dict[str, int]) -> str: +def decode_ids( + path: str, mapping: dict[str, int], *, overwrite: bool = False +) -> str: """Decode numeric IDs back to original molecule IDs. For an SDF file with numeric IDs, restores the original IDs using @@ -139,10 +147,21 @@ def decode_ids(path: str, mapping: dict[str, int]) -> str: path: Path to the input SDF file with encoded (numeric) IDs. mapping: Dictionary mapping original IDs to their numeric indices (as returned by encode_ids). + overwrite: Allow replacing an existing file at the decoded path. + Keyword-only, defaults to False: ``_out.`` is a name + Auto3D derives from `path`, so an existing file there belongs to + someone else, and ``Chem.SDWriter`` truncates on open. + ``WorkflowOrchestrator`` calls this on the combined output inside a + job directory it created with a bare ``mkdir()``, where the derived + name is always free. Returns: Path to the new SDF file with decoded IDs (adds '_out' suffix). + Raises: + ConfigurationError: A file already exists at the decoded path and + `overwrite` is False. + Example: >>> mapping = {'mol_A': 0, 'mol_B': 1} >>> output_path = decode_ids("encoded_3d.sdf", mapping) @@ -155,6 +174,8 @@ def decode_ids(path: str, mapping: dict[str, int]) -> str: stem_parts = path_obj.stem.split("_")[:-2] new_stem = "_".join(stem_parts) + "_out" new_path = path_obj.parent / f"{new_stem}.{extension}" + # Before the supplier is opened, so a refusal costs nothing. + check_output_overwrite(new_path, overwrite) suppl = Chem.SDMolSupplier(path, removeHs=False) with Chem.SDWriter(str(new_path)) as w: diff --git a/src/Auto3D/ranking.py b/src/Auto3D/ranking.py index 694513f0..b10086d9 100644 --- a/src/Auto3D/ranking.py +++ b/src/Auto3D/ranking.py @@ -6,8 +6,9 @@ from rdkit import Chem from Auto3D.config import SELECTOR_FIELDS, check_selectors_mutually_exclusive +from Auto3D.constants import DEFAULT_ENERGY_CLUSTER_WINDOW from Auto3D.exceptions import ConfigurationError, InputValidationError -from Auto3D.filtering import filter_unique, filter_unique_optimized +from Auto3D.filtering import FilterResult, filter_conformers from Auto3D.utils.connectivity import check_connectivity from Auto3D.utils.convergence import converged_or_unfiltered, has_convergence_flag from Auto3D.utils.energy import ( @@ -59,6 +60,59 @@ def species_id(name: str) -> str: return name.strip().rsplit("_", 2)[0].strip() +#: Selector field (from :data:`Auto3D.config.SELECTOR_FIELDS`) -> the +#: ``ConformerRanker`` method that implements it. +#: +#: ``run`` used to dispatch with a hand-written ``if self.k: ... elif +#: self.window: ...``, which meant a third selector added to ``SELECTOR_FIELDS`` +#: was accepted by ``Auto3DOptions``, accepted by ``CLIConfig``, accepted by +#: ``check_selectors_mutually_exclusive`` -- and then silently ignored here, +#: falling through to the "Parameter k or window needs to be specified" error +#: even though the user had specified one. The registry plus the import-time +#: equality check below turns that into an ImportError at the moment the field +#: list and this table disagree. +_SELECTORS: dict[str, str] = {"k": "top_k", "window": "top_window"} + + +def _verify_selector_registry( + registry: dict[str, str], fields: tuple[str, ...], cls: type +) -> None: + """Raise unless ``registry`` covers ``fields`` with methods that exist. + + Called at import time, and a plain function rather than inline module code + so the failure it guards against can be tested by calling it with a bad + registry -- instead of reloading the module and hoping the reload's own + machinery does not mask the error. + + Args: + registry: The selector-field -> method-name table (:data:`_SELECTORS`). + fields: The authoritative field list (``config.SELECTOR_FIELDS``). + cls: The class expected to implement each mapped method. + + Raises: + ImportError: The registry's field set differs from ``fields``, or a + mapped method is missing from ``cls``. + """ + if set(registry) != set(fields): + raise ImportError( + "Auto3D.ranking._SELECTORS is out of step with " + f"Auto3D.config.SELECTOR_FIELDS: registry has {sorted(registry)}, " + f"config declares {sorted(fields)}. config.py owns the field list; " + "add the missing selector here and give ConformerRanker a method " + "that implements it (a selector nothing dispatches on is accepted " + "by every config class and then silently ignored)." + ) + for field, method in registry.items(): + # Checked as well as the names, because a typo in a method name passes + # the set comparison above and then raises AttributeError from inside + # `run` -- after the whole input has been read and grouped. + if not callable(getattr(cls, method, None)): + raise ImportError( + f"Auto3D.ranking._SELECTORS maps selector {field!r} to " + f"{method!r}, which is not a method of {cls.__name__}." + ) + + class ConformerRanker: """Select conformers that satisfy user-defined energy criteria. @@ -72,12 +126,12 @@ class ConformerRanker: k: Select top-k structures per SMILES. False to disable. window: Select structures within window kcal/mol of minimum. False to disable. - use_optimized_filtering: If True (default), use energy-clustered - RMSD filtering which has O(n*k) complexity instead of O(n^2). - Set to False to use legacy O(n^2) filtering. - energy_cluster_window: Energy window in eV for clustering molecules - before RMSD comparison. Only used when use_optimized_filtering - is True. Default is 0.1 eV. + energy_cluster_window: Energy width in eV below which molecules are + RMSD-compared to each other. A performance knob only -- see + ``Auto3D.filtering.filter_conformers``, which is the single + conformer filter now that ``use_optimized_filtering`` and the + legacy all-pairs implementation it selected are gone. Defaults to + ``Auto3D.constants.DEFAULT_ENERGY_CLUSTER_WINDOW``. overwrite: Allow writing over an existing `out_path`. Defaults to True, which is the historical behavior every caller was written against -- including Auto3D's own pipeline, which always writes @@ -92,8 +146,7 @@ def __init__( threshold: float, k: int | bool = False, window: float | bool = False, - use_optimized_filtering: bool = True, - energy_cluster_window: float = 0.1, + energy_cluster_window: float = DEFAULT_ENERGY_CLUSTER_WINDOW, overwrite: bool = True, ) -> None: # Same C14 guard the three API entry points run. ConformerRanker is a @@ -117,29 +170,62 @@ def __init__( self.threshold = threshold self.k = k self.window = window - self.use_optimized_filtering = use_optimized_filtering self.energy_cluster_window = energy_cluster_window - - def _filter_mols(self, mols: list[Chem.Mol]) -> list[Chem.Mol]: + #: Run-level drop tally, keyed by ``Auto3D.filtering.DROP_REASONS``. + #: ``top_k``/``top_window`` add their per-species counts here (see + #: ``_account``) so ``run``'s "selected 0 structures" warning can state + #: the actual reasons instead of listing the ones it might have been. + #: Reset at the top of ``run``; harmless (and still accurate) for a + #: caller who invokes ``top_k``/``top_window`` directly. + self._drop_totals: dict[str, int] = {} + + def _account(self, result: FilterResult) -> FilterResult: + """Add ``result``'s drop counts to this run's tally, and return it.""" + for reason, count in result.dropped.items(): + if count: + self._drop_totals[reason] = self._drop_totals.get(reason, 0) + count + return result + + def _filter_mols(self, mols: list[Chem.Mol]) -> FilterResult: """Filter molecules to remove duplicates based on RMSD. - Uses either the optimized energy-clustered approach or the legacy - O(n^2) comparison, depending on the use_optimized_filtering setting. - Args: mols: List of RDKit Mol objects to filter. Returns: - List of unique molecules, sorted by energy. + A ``FilterResult``: unique molecules sorted by energy, plus the + per-reason drop counts the "nothing survived" messages below need. """ - if self.use_optimized_filtering: - return filter_unique_optimized( - mols, - rmsd_threshold=self.threshold, - energy_cluster_window=self.energy_cluster_window, - ) + return filter_conformers( + mols, + rmsd_threshold=self.threshold, + energy_cluster_window=self.energy_cluster_window, + ) + + def _log_nothing_selected(self, name: str, result: FilterResult) -> None: + """Say WHY a species contributed no structure to the output. + + This used to be an unconditional "No structure converged for X." -- + which is what a reader saw when every conformer of X was dropped for + *stereochemistry* (an optimization that inverted a center, so the + geometry no longer matches the title) or for *connectivity* (a + structure that fell apart). Both point at the input or the chemistry, + and both were reported as an optimizer convergence problem, sending the + reader to ``--opt-steps`` and ``--convergence-threshold`` for something + neither could fix. + + The literal wording survives for the one case where it was true -- + convergence being the sole reason -- because it is the message users + and their log-scraping scripts have matched on since Auto3D 1.x. When + the sole reason is convergence there is nothing to add; when it is not, + the reasons are named. + """ + if set(result.dropped) <= {"unconverged"}: + logger.info(f"No structure converged for {name}.") else: - return filter_unique(mols, crit=self.threshold) + logger.info( + "No structure selected for %s: %s.", name, result.summary() + ) def top_k(self, df_group: pd.DataFrame, k: int = 1) -> list[Chem.Mol]: """Select top-k lowest-energy structures from a group. @@ -159,28 +245,43 @@ def top_k(self, df_group: pd.DataFrame, k: int = 1) -> list[Chem.Mol]: df2 = df_group.sort_values(by=['energies']) # Optimization: when k=1, skip the expensive RMSD dedup but still - # apply connectivity validation. Return the lowest-energy conformer - # that passes check_connectivity (no broken/formed bonds); if none - # pass, return an empty list. + # apply the validity checks. Return the lowest-energy conformer that + # passes them; if none pass, return an empty list. + # + # The predicates and the drop-reason names are deliberately the same + # ones `filter_conformers` uses, so the diagnostic a user gets for an + # empty selection does not depend on which k they asked for. (Dedup is + # the only thing skipped -- it cannot empty a non-empty set anyway.) + # The loop stops at the first survivor, so the counts are partial when + # something IS selected and complete in exactly the case that produces + # a message. if k == 1: out_mols = [] + dropped: dict[str, int] = {} for mol in df2["mols"]: - if stereo_preserved(mol) and check_connectivity(mol): + if not converged_or_unfiltered(mol): + reason = "unconverged" + elif not stereo_preserved(mol): + reason = "stereochemistry" + elif not check_connectivity(mol): + reason = "connectivity" + else: out_mols = [mol] break + dropped[reason] = dropped.get(reason, 0) + 1 + result = self._account(FilterResult(kept=out_mols, dropped=dropped)) else: - out_mols_ = self._filter_mols(list(df2["mols"])) - if k < len(out_mols_): - out_mols = out_mols_[:k] - else: - out_mols = out_mols_ + result = self._account(self._filter_mols(list(df2["mols"]))) + # Truncation to k is selection, not a filter drop: those conformers + # are valid and unique, they just lost the ranking, so they are not + # counted among `dropped` (nothing is "missing" to explain). + out_mols = result.kept[:k] if k < len(result.kept) else result.kept if len(out_mols) == 0: # names[0] is already the group's species id (see species_id() # above) -- no further splitting here, or a disambiguated id # like "KEY_2" would misreport as "KEY" in this message. - name = names[0].strip() - logger.info(f"No structure converged for {name}.") + self._log_nothing_selected(names[0].strip(), result) else: #Adding relative energies # E_tot is stored in Hartree; E_rel(eV) is, as its name says, eV. @@ -210,21 +311,20 @@ def top_window(self, df_group: pd.DataFrame, window: float = 1.0) -> list[Chem.M raise ValueError(f"All molecules must have the same name, got: {set(names)}") df2 = df_group.sort_values(by=['energies']) - out_mols_ = self._filter_mols(list(df2['mols'])) + result = self._account(self._filter_mols(list(df2['mols']))) out_mols = [] - if len(out_mols_) == 0: + if len(result.kept) == 0: # names[0] is already the group's species id -- see the note in # top_k above. - name = names[0].strip() - logger.info(f"No structure converged for {name}.") + self._log_nothing_selected(names[0].strip(), result) else: # `window` was converted from kcal/mol to eV above, and E_tot is # stored in Hartree, so both sides of the comparison are eV here. # Reading the Hartree number as if it were eV is what made the # window 27.2x too wide for an opt_geometry-produced input. - ref_energy = e_tot_ev(out_mols_[0]) - for mol in out_mols_: + ref_energy = e_tot_ev(result.kept[0]) + for mol in result.kept: my_energy = e_tot_ev(mol) rel_energy = my_energy - ref_energy if rel_energy <= window: @@ -232,6 +332,16 @@ def top_window(self, df_group: pd.DataFrame, window: float = 1.0) -> list[Chem.M out_mols.append(mol) else: break + # The window is the one drop reason this method owns, and it is + # merged into the same tally as the filter's own counts rather than + # reported separately, so `run`'s summary is one accounting of the + # whole selection. `break` above is safe because `kept` ascends in + # energy, so every remaining conformer is outside the window too. + n_outside = len(result.kept) - len(out_mols) + if n_outside: + self._account( + FilterResult(kept=out_mols, dropped={"energy_window": n_outside}) + ) return out_mols def run(self) -> list[Chem.Mol]: @@ -271,9 +381,13 @@ def run(self) -> list[Chem.Mol]: {name: getattr(self, name, None) for name in SELECTOR_FIELDS} ) results = [] + # Fresh tally per run, so a ranker reused for a second file does not + # report the first file's drops (see _account). + self._drop_totals = {} mols, names, energies = [], [], [] n_records = 0 + n_unparsed = 0 n_unconverged = 0 n_unflagged = 0 # Context-managed so the SDF file handle is released promptly rather than @@ -281,6 +395,7 @@ def run(self) -> list[Chem.Mol]: with Chem.SDMolSupplier(self.input_path, removeHs=False) as supplier: for position, mol in enumerate(supplier): if mol is None: + n_unparsed += 1 logger.warning( "Skipping record %d of %s: RDKit could not parse it.", position, self.input_path, @@ -327,10 +442,18 @@ def run(self) -> list[Chem.Mol]: for group_name in groups.indices: group = groups.get_group(group_name) - if self.k: - top_results = self.top_k(group, self.k) - elif self.window: - top_results = self.top_window(group, self.window) + # Dispatch through the registry rather than a hand-written + # if/elif chain on the field names, so the set of selectors this + # method honors is the set `Auto3D.config` declares -- checked for + # equality at import (see _SELECTORS). Iterating SELECTOR_FIELDS + # rather than the dict also means config.py owns the PRECEDENCE, + # not this module's dict-literal order; the mutual-exclusivity + # check above has already refused more than one anyway. + for field in SELECTOR_FIELDS: + value = getattr(self, field, None) + if value: + top_results = getattr(self, _SELECTORS[field])(group, value) + break else: raise ConfigurationError('Parameter k or window needs to be ' 'specified. Append "--k=1" if you' @@ -342,12 +465,28 @@ def run(self) -> list[Chem.Mol]: # returns []. WARNING, not INFO: `logging.lastResort` puts WARNING # and above on stderr even for a caller who never ran # configure_logging, which is every direct API caller. + # + # The reasons are the ones that actually fired. This used to read + # "N record(s) are marked Converged=false and the rest were dropped + # by the connectivity, stereochemistry or energy-window filters" -- + # a hand-maintained list of everything it MIGHT have been, which + # left the reader to work out which, and which had to be edited by + # hand every time a filter was added. `_drop_totals` is the tally + # `top_k`/`top_window` fed as they ran; the two counts below are + # the drops `run` itself made, before grouping, which never reach a + # filter and so are not double-counted. + totals = dict(self._drop_totals) + for reason, count in ( + ("unparsed", n_unparsed), ("unconverged", n_unconverged), + ): + if count: + totals[reason] = totals.get(reason, 0) + count + summary = FilterResult(kept=[], dropped=totals).summary() logger.warning( "Selected 0 structures from %d record(s) in %s, so %s is " - "empty: %d record(s) are marked Converged=false and the rest " - "were dropped by the connectivity, stereochemistry or " - "energy-window filters.", - n_records, self.input_path, self.out_path, n_unconverged, + "empty: %s.", + n_records, self.input_path, self.out_path, + summary or "no filter reported a reason", ) with Chem.SDWriter(self.out_path) as f: @@ -372,5 +511,11 @@ def run(self) -> list[Chem.Mol]: return results +# Runs at import: a selector declared in config.py with nothing wired to it +# here, or wired to a method that does not exist, fails the import rather than +# being discovered by a user whose `--k=1` was silently ignored. +_verify_selector_registry(_SELECTORS, SELECTOR_FIELDS, ConformerRanker) + + # Backward compatibility alias ranking = ConformerRanker diff --git a/src/Auto3D/tautomer.py b/src/Auto3D/tautomer.py index 438941e1..a38add99 100644 --- a/src/Auto3D/tautomer.py +++ b/src/Auto3D/tautomer.py @@ -10,13 +10,20 @@ from Auto3D.exceptions import ConfigurationError from Auto3D.utils.energy import e_tot_hartree, hartree2kcalpermol from Auto3D.utils.logging_config import get_logger +from Auto3D.utils.output_guard import check_output_overwrite logger = get_logger(__name__) __all__ = ["select_tautomers", "get_stable_tautomers"] -def select_tautomers(sdf: str, k: int | None = None, window: float | None = None) -> str: +def select_tautomers( + sdf: str, + k: int | None = None, + window: float | None = None, + *, + overwrite: bool = False, +) -> str: """Select and Write the top-k or E <= window tautomers for each input SMILES Only k or window needs to be specified, NOT both. @@ -24,18 +31,23 @@ def select_tautomers(sdf: str, k: int | None = None, window: float | None = None Output: the path of the low-energy tautomer 3D conformers - Warning: - This function writes ``/_top_tautomers.sdf`` with no - overwrite gate: ``Chem.SDWriter`` truncates on open, so a direct call - such as ``select_tautomers("/data/results.sdf", k=1)`` replaces any - existing ``/data/results_top_tautomers.sdf``. Every route Auto3D - itself takes is safe -- ``get_stable_tautomers`` passes ``main()``'s - output, which lives in a job directory created fresh for that run, and - ``auto3d tautomers`` additionally gates its ``-o`` with - ``check_output_overwrite`` -- so this is a hazard for direct API - callers only, the same residual class as - ``Auto3D.utils.smi_io.smiles2smi`` and ``Auto3D.id_mapping.decode_ids``. See - ``docs/superpowers/follow-ups-after-4.0.0-remediation.md``. + Args: + sdf: Path to the SDF this function reads, normally ``main()``'s output. + k: Keep the top-k tautomers per input molecule. + window: Keep tautomers within this kcal/mol window of the most stable + one. Mutually exclusive with `k`. + overwrite: Allow replacing an existing + ``/_top_tautomers.sdf``. Keyword-only, and + defaults to False because that name is one Auto3D *derives* rather + than one the caller chose: ``Chem.SDWriter`` truncates on open, so + ``select_tautomers("/data/results.sdf", k=1)`` used to replace an + existing ``/data/results_top_tautomers.sdf`` with this call's + selection, silently. Auto3D's own routes are unaffected -- + ``get_stable_tautomers`` passes ``main()``'s output, which lives in + a job directory created fresh for that run, and ``auto3d + tautomers`` additionally gates its ``-o`` with + ``check_output_overwrite`` -- so the default only ever fires for a + direct API caller pointing at an occupied path. Note: Tautomers are ranked by the optimized NNP *electronic* energy (``E_tot``) @@ -48,7 +60,8 @@ def select_tautomers(sdf: str, k: int | None = None, window: float | None = None Raises: ConfigurationError: If both k and window are given, if neither is - given, or if k < 1. ``auto3d tautomers`` already rejects the + given, if k < 1, or if the derived output path exists and + `overwrite` is False. ``auto3d tautomers`` already rejects the both-given case in ``execute_tautomers`` before calling this function, but ``select_tautomers``/``get_stable_tautomers`` are also public Python API entry points that can be called directly, @@ -60,6 +73,16 @@ def select_tautomers(sdf: str, k: int | None = None, window: float | None = None if (k is not None) and (k < 1): raise ConfigurationError(f"tauto_k must be >= 1, got {k}") + # Resolved and gated up front, before the input is read and grouped: + # refusing after all the work is done costs the user the run for nothing, + # and the writer at the bottom truncates the moment it opens. + folder = os.path.dirname(sdf) + # splitext (not split(".")) so an input like 'mol.v2.sdf' keeps 'mol.v2' + # instead of collapsing to 'mol' and risking output collisions. + stem = os.path.splitext(os.path.basename(sdf))[0].strip() + output_path = os.path.join(folder, stem + "_top_tautomers.sdf") + check_output_overwrite(output_path, overwrite) + supplier = Chem.SDMolSupplier(sdf, removeHs=False) mols = [m for m in supplier if m is not None] for mol in mols: @@ -101,14 +124,7 @@ def select_tautomers(sdf: str, k: int | None = None, window: float | None = None else: raise ConfigurationError("Either k OR window needs to be specified") results += out_mols - - folder = os.path.dirname(sdf) - # splitext (not split(".")) so an input like 'mol.v2.sdf' keeps 'mol.v2' - # instead of collapsing to 'mol' and risking output collisions. - stem = os.path.splitext(os.path.basename(sdf))[0].strip() - basename = stem + "_top_tautomers.sdf" - output_path = os.path.join(folder, basename) with Chem.SDWriter(output_path) as w: for mol in results: w.write(mol) diff --git a/src/Auto3D/utils/convergence.py b/src/Auto3D/utils/convergence.py index 68b2505b..b8570922 100644 --- a/src/Auto3D/utils/convergence.py +++ b/src/Auto3D/utils/convergence.py @@ -4,8 +4,7 @@ ``Converged`` is written by :mod:`Auto3D.batch_opt.batchopt` for every record it optimizes -- ``"True"`` or ``"False"`` -- and read by the three filters that decide which conformers survive (:class:`Auto3D.ranking.ConformerRanker`, -:func:`Auto3D.filtering.filter_unique_optimized`, -:func:`Auto3D.filtering.filter_unique`). +:func:`Auto3D.filtering.filter_conformers`). All three used to read it as:: diff --git a/src/Auto3D/utils/energy.py b/src/Auto3D/utils/energy.py index 76fc2787..3ebaecfc 100644 --- a/src/Auto3D/utils/energy.py +++ b/src/Auto3D/utils/energy.py @@ -10,8 +10,7 @@ ``E_tot`` in eV and ``ASE/geometry.opt_geometry`` converted the same tag to Hartree afterwards, so the identical property name carried two units depending on which entry point produced the file -- and the five in-package consumers -(``ranking``, ``filtering.filter_unique_optimized``, -``filtering.filter_unique``) all hard-coded +(``ranking``, ``filtering.filter_conformers``) all hard-coded eV. Feeding an ``opt_geometry`` output straight to ``ConformerRanker(window=2.0)`` therefore opened a window 27.2x too wide, kept 3 conformers where 2 belong, reported ``E_rel`` 0.037 kcal/mol where the truth diff --git a/src/Auto3D/utils/geometry.py b/src/Auto3D/utils/geometry.py index ca8c5018..2f89ae9f 100644 --- a/src/Auto3D/utils/geometry.py +++ b/src/Auto3D/utils/geometry.py @@ -82,7 +82,7 @@ def get_rmsd(mol1: Chem.Mol, mol2: Chem.Mol, remove_hs: bool = True) -> float: The RMSD value in Angstroms. Returns ``float("inf")`` if alignment fails (e.g., due to atom mismatch). An incomparable pair is treated as "distinct" rather than "identical", which is the same convention used - by ``filter_unique``; a downstream ``rmsd < threshold`` check therefore + by the duplicate filter; a downstream ``rmsd < threshold`` check therefore keeps the structure instead of dropping it as a false duplicate. Example: @@ -106,6 +106,6 @@ def get_rmsd(mol1: Chem.Mol, mol2: Chem.Mol, remove_hs: bool = True) -> float: # Temporary bug fix for https://github.com/rdkit/rdkit/issues/6826 rmsd = rdMolAlign.GetBestRMS(mol1_proc, mol2_proc) except RuntimeError: - # Incomparable pair: treat as distinct (inf), matching filter_unique. + # Incomparable pair: treat as distinct (inf), as the filter requires. rmsd = float("inf") return float(rmsd) diff --git a/src/Auto3D/utils/output_guard.py b/src/Auto3D/utils/output_guard.py index 0760452e..4d778ecb 100644 --- a/src/Auto3D/utils/output_guard.py +++ b/src/Auto3D/utils/output_guard.py @@ -8,8 +8,9 @@ ``Auto3D.exceptions``, which is what lets the ``.smi``/``.sdf`` writers under ``utils/`` and the top-level ID/layout helpers gate their own output. -``utils/validation.py`` still re-exports both names, because ``SPE``, -``ASE.thermo`` and ``cli.commands.properties`` import them from there. +This is the only path to these two names. ``utils/validation.py`` re-exported +them for one release while call sites moved over; that re-export is gone, so +importing them from ``validation`` now fails rather than quietly working. """ from __future__ import annotations diff --git a/src/Auto3D/utils/smi_io.py b/src/Auto3D/utils/smi_io.py index 4781705d..38b32523 100644 --- a/src/Auto3D/utils/smi_io.py +++ b/src/Auto3D/utils/smi_io.py @@ -23,6 +23,7 @@ from Auto3D.exceptions import InputValidationError from Auto3D.utils.logging_config import get_logger +from Auto3D.utils.output_guard import check_output_overwrite logger = get_logger(__name__) @@ -86,7 +87,7 @@ def iter_smi_records(path, *, on_malformed="skip"): yield line_no, parts[0], parts[1] -def smiles2smi(smiles: list[str], path: str) -> str: +def smiles2smi(smiles: list[str], path: str, *, overwrite: bool = True) -> str: """Convert a list of SMILES strings to a .smi file with InChIKey IDs. Each SMILES string is converted to a molecule, and its InChIKey is computed @@ -96,10 +97,22 @@ def smiles2smi(smiles: list[str], path: str) -> str: Args: smiles: List of SMILES strings to convert. path: Output file path for the .smi file. + overwrite: Allow replacing an existing file at `path`. Keyword-only. + Defaults to **True**, unlike the writers that derive their own + output name (``id_mapping.encode_ids``/``decode_ids``, + ``tautomer.select_tautomers``): the caller named this file, so + refusing it would gate a decision the caller already made. + ``smiles2mols`` writes into a fresh ``TemporaryDirectory`` on every + call and relies on that. Pass False to make the write refuse an + occupied path. Returns: The output file path. + Raises: + ConfigurationError: `path` exists and `overwrite` is False. + InputValidationError: A SMILES string could not be parsed by RDKit. + Example: >>> smiles2smi(["CCO", "CCC"], "molecules.smi") 'molecules.smi' @@ -107,6 +120,11 @@ def smiles2smi(smiles: list[str], path: str) -> str: # CCO LFQSCWFLJHTTHZ-UHFFFAOYSA-N # CCC ATUOYWHBWRKTHZ-UHFFFAOYSA-N """ + # Checked before any SMILES is parsed: refusing after the work is done is a + # worse experience for no gain, and `open(path, "w+")` below truncates on + # open, so by the time the write starts the file is already gone. + check_output_overwrite(path, overwrite) + lines = [] seen_ids: dict[str, int] = {} for idx, smi in enumerate(smiles): diff --git a/src/Auto3D/utils/stereochemistry.py b/src/Auto3D/utils/stereochemistry.py index f6d7e2af..14c2ba51 100644 --- a/src/Auto3D/utils/stereochemistry.py +++ b/src/Auto3D/utils/stereochemistry.py @@ -18,6 +18,7 @@ from Auto3D.utils.atomic_io import atomic_write_path from Auto3D.utils.logging_config import get_logger +from Auto3D.utils.smi_io import iter_smi_records logger = get_logger(__name__) @@ -252,12 +253,16 @@ def remove_enantiomers(inpath: str, out: str) -> dict[str, list[str]]: If enantiomer removal fails for a molecule, all original SMILES are kept and a message is printed. """ - with open(inpath) as f: - data = f.readlines() - + # iter_smi_records is the single parser for this format (M59). The loop + # this replaced built `vals = line.split()` and indexed vals[0]/vals[1] + # with no guard, so a blank line in `inpath` raised a bare IndexError; + # "raise" keeps that fail-fast behavior (as an InputValidationError naming + # the line, not a crash) for any genuinely malformed row while blank lines + # -- which should never appear in a file this module itself wrote, but + # previously crashed the whole function if one did -- are now skipped like + # every other consumer of this format. smiles: dict[str, list[str]] = defaultdict(lambda: []) - for line in data: - vals = line.split() + for _line_no, smi, mol_id in iter_smi_records(inpath, on_malformed="raise"): # Strip only the trailing isomer-index component write_enumerated_smi # appends (rsplit, maxsplit=1), not everything after the first # underscore: an id like "KEY_2" -- smiles2smi's disambiguation of a @@ -265,7 +270,7 @@ def remove_enantiomers(inpath: str, out: str) -> dict[str, list[str]]: # specifically so it is not dropped -- must survive this grouping # intact, or it silently merges back onto "KEY" here before ranking # ever sees it (M17). - smi, name = vals[0].strip(), vals[1].strip().rsplit("_", 1)[0].strip() + name = mol_id.rsplit("_", 1)[0] smiles[name].append(smi) for key, values in smiles.items(): @@ -407,47 +412,35 @@ def create_enantiomer(smi: str) -> str: 'C[C@@H](O)F' """ stereo_info = get_stereo_info(smi) - new_smi = "" keys = list(stereo_info.keys()) - if len(keys) == 0: + if not keys: # No stereo centers to invert return smi - if len(keys) == 1: - key = keys[0] + + # Single pass with a cursor, one key at a time: emit smi[cursor:key] + # verbatim, then the inverted marker, then advance the cursor past the + # marker just consumed. The tail after the last key is emitted once, + # after the loop, from the final cursor position. + # + # This replaces two copies of the same logic that used to disagree in + # shape: a `len(keys) == 1` branch handled one stereo center, and a + # general loop handled two or more by reading `key2`/`val2` (set inside + # the loop's `else` branch) *after* the `for` -- correct only because + # Python does not scope loop variables to the loop body, so `key2` still + # held the last iteration's value. That accident breaks the moment + # `keys` can be empty going into the tail read, which is exactly what the + # `len(keys) == 1` branch existed to avoid (M60). + inverted = {"@": "@@", "@@": "@"} + new_smi = "" + cursor = 0 + for key in keys: val = stereo_info[key] - if val == "@": - new_smi += smi[:key] - new_smi += "@@" - new_smi += smi[(key + 1) :] - elif val == "@@": - new_smi += smi[:key] - new_smi += "@" - new_smi += smi[(key + 2) :] - else: + if val not in inverted: raise ValueError("Invalid %s" % smi) - return new_smi - - for i in range(len(keys)): - if i == 0: - key = keys[i] - new_smi += smi[:key] - else: - key1 = keys[i - 1] - key2 = keys[i] - val1 = stereo_info[key1] - if val1 == "@": - new_smi += "@@" - new_smi += smi[int(key1 + 1) : key2] - elif val1 == "@@": - new_smi += "@" - new_smi += smi[int(key1 + 2) : key2] - val2 = stereo_info[key2] - if val2 == "@": - new_smi += "@@" - new_smi += smi[int(key2 + 1) :] - elif val2 == "@@": - new_smi += "@" - new_smi += smi[int(key2 + 2) :] + new_smi += smi[cursor:key] + new_smi += inverted[val] + cursor = key + len(val) + new_smi += smi[cursor:] return new_smi @@ -499,16 +492,20 @@ def amend_configuration(smis: str) -> dict[str, list[str]]: For input "N=C1OC(CN2CC(C)OC(C)C2)CN1", if some stereo configurations are missing, this function will attempt to add them. """ - with open(smis) as f: - data = f.readlines() + # iter_smi_records is the single parser for this format (M59). The + # `tuple(line.strip().split())` this replaced required EXACTLY 2 tokens + # (a third raised "too many values to unpack") and had no blank/comment + # handling at all; "raise" preserves the fail-fast behavior for a + # genuinely malformed row (now a named InputValidationError instead of a + # bare ValueError) while tolerating a trailing extra column, matching + # every other consumer of this format. dct: dict[str, list[str]] = defaultdict(lambda: []) - for line in data: - smi, idx = tuple(line.strip().split()) + for _line_no, smi, mol_id in iter_smi_records(smis, on_malformed="raise"): # See the matching note in remove_enantiomers: strip only the # trailing isomer-index component, not everything after the first # underscore, so a disambiguated id like "KEY_2" is not merged back # onto "KEY" here (M17). - idx = idx.rsplit("_", 1)[0].strip() + idx = mol_id.rsplit("_", 1)[0] dct[idx].append(smi) for key in dct.keys(): diff --git a/src/Auto3D/utils/validation.py b/src/Auto3D/utils/validation.py index af2656a9..deed6858 100644 --- a/src/Auto3D/utils/validation.py +++ b/src/Auto3D/utils/validation.py @@ -22,16 +22,7 @@ ModelLoadError, ) from Auto3D.utils.logging_config import get_logger - -# The two output guards moved to the leaf module ``utils/output_guard.py`` so -# that a writer needing only them does not import this module's torch. They are -# re-exported here because ``SPE``, ``ASE.thermo`` and -# ``cli.commands.properties`` import them from this path; new call sites should -# name ``Auto3D.utils.output_guard`` directly. -from Auto3D.utils.output_guard import ( # noqa: F401 - check_output_not_input, - check_output_overwrite, -) +from Auto3D.utils.smi_io import iter_smi_records from Auto3D.utils.stereochemistry import count_unspecified_stereo if TYPE_CHECKING: @@ -291,31 +282,15 @@ def check_smi_format(args: Any) -> tuple[bool, list[str]]: """ ANI = True - smiles_all = [] - with open(args.path) as f: - data = f.readlines() - for line in data: - if line.isspace(): - continue - # Skip '#'-prefixed comment lines, matching cli.commands.validate. - # validate_smiles_file and file_ops.iter_smi_records -- all three must - # agree on what counts as a comment vs. data, or `auto3d validate` - # would approve a file this function then rejects (M25). - if line.lstrip().startswith("#"): - continue - # Tolerate ragged rows the way the rest of the pipeline does: the chunk - # loader reads only the first two whitespace columns (usecols=[0, 1]), so - # trailing tokens (e.g. an inline comment column) must not be rejected - # here. split() never yields empty tokens, so a present SMILES/ID is - # guaranteed non-empty. - parts = line.split() - if len(parts) < 2: - raise InputValidationError( - "Each non-blank line must contain a SMILES and an ID separated by " - f"whitespace, but got: {line.strip()!r}" - ) - smiles = parts[0] # parts[1] is the ID; its presence is enforced above - smiles_all.append(smiles) + # iter_smi_records is the single parser for this format (M59): it already + # skips blank lines and '#'-prefixed comments the same way + # cli.commands.validate.validate_smiles_file does, and its on_malformed + # ("raise") gives the same InputValidationError this loop used to raise by + # hand for a line missing an ID, so `auto3d validate` and this check + # cannot silently disagree about what a well-formed line looks like + # (M25). The parser also tolerates ragged rows (extra whitespace columns + # beyond SMILES+ID), matching the chunk loader's usecols=[0, 1]. + smiles_all = [smiles for _line_no, smiles, _id in iter_smi_records(args.path, on_malformed="raise")] logger.info(f"\tThere are {len(smiles_all)} SMILES in the input file {args.path}.") logger.info("\tAll SMILES and IDs are valid.") diff --git a/tests/test_chunk_manager.py b/tests/test_chunk_manager.py index e5a0c960..fe83549f 100644 --- a/tests/test_chunk_manager.py +++ b/tests/test_chunk_manager.py @@ -395,6 +395,147 @@ def test_chunk_size_clamped_to_at_least_one(self, tmp_path): assert 1 <= len(chunk_info) <= 4 +class TestSmiParserCrossAgreement: + """M59: chunk_manager's pandas-based `.smi` reader stays separate from + ``iter_smi_records`` for performance (avoiding a pure-Python per-line loop + over what can be a very large file), but the two must not silently drift + apart on the input they both exist to read: an already-ID-encoded + intermediate `.smi` file with no blank lines and no comments (chunk_manager + calls this "encode_ids semantics" in prepare_chunks' docstring). If either + parser's tokenizing/whitespace-splitting rule changes, this test should + catch the divergence. + """ + + def test_prepare_chunks_agrees_with_iter_smi_records_on_well_formed_input( + self, tmp_path + ): + """Both parsers must extract the same (smiles, id) pairs, in the same + order, from a well-formed encoded .smi file.""" + from Auto3D.utils.smi_io import iter_smi_records + + rows = [ + ("CCO", "0"), + ("CCCO", "1"), + ("c1ccccc1", "2"), + ("C[C@H](O)F", "3"), + ("CCN", "4"), + ] + input_file = tmp_path / "test_encoded.smi" + input_file.write_text("".join(f"{smi} {mol_id}\n" for smi, mol_id in rows)) + + # memory * capacity is chosen well above len(rows) so chunk_size never + # forces a split: exactly one chunk, holding every row in original + # order, so the chunk file's content is directly comparable to + # iter_smi_records' output order (round-robin distribution across + # multiple chunks would otherwise reorder rows relative to the input). + config = Auto3DOptions( + path=str(tmp_path / "test.smi"), k=1, memory=1, capacity=1000 + ) + manager = ChunkManager( + config=config, + input_path=input_file, + input_format="smi", + job_dir=tmp_path, + workflow_logger=None, + ) + + chunk_info = manager.prepare_chunks() + assert len(chunk_info) == 1, "test assumes no chunk splitting" + + chunk_path, _ = chunk_info[0] + pandas_rows = [ + tuple(line.split()) + for line in Path(chunk_path).read_text().splitlines() + if line.strip() + ] + iter_rows = [ + (smi, mol_id) for _line_no, smi, mol_id in iter_smi_records(str(input_file)) + ] + + assert pandas_rows == rows + assert iter_rows == rows + assert pandas_rows == iter_rows + + def test_ragged_extra_column_agreement(self, tmp_path): + """A trailing whitespace-separated column beyond SMILES+ID must be + dropped identically by both readers (chunk_manager's usecols=[0, 1] + vs. iter_smi_records taking only parts[0]/parts[1]).""" + from Auto3D.utils.smi_io import iter_smi_records + + input_file = tmp_path / "test_encoded.smi" + input_file.write_text("CCO 0 inline_comment_column\nCCCO 1\n") + + config = Auto3DOptions( + path=str(tmp_path / "test.smi"), k=1, memory=1, capacity=1000 + ) + manager = ChunkManager( + config=config, + input_path=input_file, + input_format="smi", + job_dir=tmp_path, + workflow_logger=None, + ) + chunk_info = manager.prepare_chunks() + assert len(chunk_info) == 1 + + chunk_path, _ = chunk_info[0] + pandas_rows = [ + tuple(line.split()) + for line in Path(chunk_path).read_text().splitlines() + if line.strip() + ] + iter_rows = [ + (smi, mol_id) for _line_no, smi, mol_id in iter_smi_records(str(input_file)) + ] + + assert pandas_rows == iter_rows == [("CCO", "0"), ("CCCO", "1")] + + def test_documented_divergence_comment_lines(self, tmp_path): + """Documents a real, deliberate divergence rather than papering over + it: iter_smi_records treats a '#'-prefixed line as a comment and + skips it (matching cli.commands.validate, per M25); chunk_manager's + pd.read_csv has no `comment=` parameter and reads it as a data row. + + This is not a bug this task fixes (chunk_manager stays a separate, + faster reader by design) -- it is recorded here so that a future + change to either parser's comment handling is a deliberate, + visible decision rather than a silent behavior change caught only in + production. If this test starts failing because someone taught + pd.read_csv to skip '#' lines, that is progress: update the + assertion, don't just delete the test. + """ + from Auto3D.utils.smi_io import iter_smi_records + + input_file = tmp_path / "test_encoded.smi" + input_file.write_text("CCO 0\n# 1 2\nCCCO 3\n") + + config = Auto3DOptions( + path=str(tmp_path / "test.smi"), k=1, memory=1, capacity=1000 + ) + manager = ChunkManager( + config=config, + input_path=input_file, + input_format="smi", + job_dir=tmp_path, + workflow_logger=None, + ) + chunk_info = manager.prepare_chunks() + chunk_path, _ = chunk_info[0] + pandas_rows = [ + tuple(line.split()) + for line in Path(chunk_path).read_text().splitlines() + if line.strip() + ] + iter_rows = [ + (smi, mol_id) for _line_no, smi, mol_id in iter_smi_records(str(input_file)) + ] + + # pandas reads the '#' line as data; iter_smi_records skips it. + assert pandas_rows == [("CCO", "0"), ("#", "1"), ("CCCO", "3")] + assert iter_rows == [("CCO", "0"), ("CCCO", "3")] + assert pandas_rows != iter_rows + + class TestLogging: """Tests for ChunkManager logging.""" diff --git a/tests/test_cli_errors.py b/tests/test_cli_errors.py index 84b55324..d67a6a0c 100644 --- a/tests/test_cli_errors.py +++ b/tests/test_cli_errors.py @@ -288,7 +288,7 @@ def _flat(text: str) -> str: def test_overwrite_refusal_does_not_suggest_config_init(capsys, tmp_path): """Raised by the real guard, not hand-built: this pins the raise site's choice of hint together with handle_error's presentation of it.""" - from Auto3D.utils.validation import check_output_overwrite + from Auto3D.utils.output_guard import check_output_overwrite existing = tmp_path / "precious.sdf" existing.write_text("x") diff --git a/tests/test_cli_property_commands.py b/tests/test_cli_property_commands.py index 66a58d22..0ed86409 100644 --- a/tests/test_cli_property_commands.py +++ b/tests/test_cli_property_commands.py @@ -481,7 +481,7 @@ def test_tautomers_refuses_output_equal_to_input(smi): # `auto3d energy junk.sdf --no-gpu -o precious.sdf` used to exit 0, print # "Wrote precious.sdf", and leave precious.sdf at 0 bytes. `config init` has # had -f/--force since it shipped; these four commands did not. The guard -# itself (`Auto3D.utils.validation.check_output_overwrite`) is exercised per +# itself (`Auto3D.utils.output_guard.check_output_overwrite`) is exercised per # API function in tests/test_durability.py; what is pinned here is the CLI # half -- that each command actually *passes* its flag down, which is the part # a refactor drops silently. diff --git a/tests/test_durability.py b/tests/test_durability.py index a8a3c818..ed8f3c39 100644 --- a/tests/test_durability.py +++ b/tests/test_durability.py @@ -510,7 +510,7 @@ def test_guard_compares_real_paths_not_strings(self, job_dir, monkeypatch): a guard that always raised would satisfy every test above. """ from Auto3D.exceptions import ConfigurationError - from Auto3D.utils.validation import check_output_not_input + from Auto3D.utils.output_guard import check_output_not_input sdf = job_dir / "mols.sdf" _write_sdf(sdf, ["m1"]) @@ -555,7 +555,7 @@ def test_a_hardlink_to_the_input_is_refused(self, job_dir): directly on this ext4 box, but it takes the identical samefile path. """ from Auto3D.exceptions import ConfigurationError - from Auto3D.utils.validation import check_output_not_input + from Auto3D.utils.output_guard import check_output_not_input sdf = job_dir / "mols.sdf" _write_sdf(sdf, ["a"]) @@ -572,7 +572,7 @@ def test_a_hardlink_to_the_input_is_refused(self, job_dir): def test_a_genuinely_different_output_is_still_allowed(self, job_dir): """Negative control: without this, a guard that always raised would satisfy every other test in this class.""" - from Auto3D.utils.validation import check_output_not_input + from Auto3D.utils.output_guard import check_output_not_input sdf = job_dir / "mols.sdf" _write_sdf(sdf, ["a"]) @@ -976,7 +976,7 @@ def test_conformer_ranker_refuses_an_existing_output(self, job_dir): def test_the_guard_permits_everything_it_should(self, job_dir): """Negative controls for the shared function itself.""" - from Auto3D.utils.validation import check_output_overwrite + from Auto3D.utils.output_guard import check_output_overwrite existing = job_dir / "there.sdf" existing.write_text("x") @@ -992,7 +992,7 @@ def test_the_guard_accepts_str_and_path_alike(self, job_dir): an `os.path.join` string. A guard that only handled one would silently no-op for the other.""" from Auto3D.exceptions import ConfigurationError - from Auto3D.utils.validation import check_output_overwrite + from Auto3D.utils.output_guard import check_output_overwrite existing = job_dir / "there.sdf" existing.write_text("x") diff --git a/tests/test_filter_unique.py b/tests/test_filter_unique.py deleted file mode 100644 index e3ce2d35..00000000 --- a/tests/test_filter_unique.py +++ /dev/null @@ -1,268 +0,0 @@ -#!/usr/bin/env python -"""Tests for the legacy all-pairs ``Auto3D.filtering.filter_unique``. - -Kept as its own file because ``filter_unique`` is scheduled for removal once -``filter_unique_optimized`` is the single filter; deleting it then means -deleting this file, not surgery on a shared one. -""" -from __future__ import annotations - -import pytest # noqa: F401 -from rdkit import Chem -from rdkit.Chem import AllChem - -from Auto3D.filtering import filter_unique - - -class TestFilterUnique: - """Test the filter_unique function for RMSD-based duplicate filtering.""" - - def test_filter_identical_conformers(self): - """Test that identical conformers are filtered to one.""" - - mol = Chem.MolFromSmiles("CCO") - mol = Chem.AddHs(mol) - AllChem.EmbedMolecule(mol, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol) - mol.SetProp("Converged", "true") - - # Create identical copies - mol2 = Chem.Mol(mol) - mol2.SetProp("Converged", "true") - - mols = [mol, mol2] - unique_mols = filter_unique(mols, crit=0.3) - - # Should only keep one - assert len(unique_mols) == 1 - - def test_same_geometry_different_energy_kept(self): - """Identical geometry but distinct E_tot must be kept (energy guard). - - Heavy-atom RMSD ~= 0 but the two are distinct minima (the O-H rotamer - case); the energy guard must stop them collapsing into one. - """ - - mol = Chem.MolFromSmiles("CCO") - mol = Chem.AddHs(mol) - AllChem.EmbedMolecule(mol, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol) - mol.SetProp("Converged", "true") - mol.SetProp("E_tot", "-10.0") - - mol2 = Chem.Mol(mol) # identical geometry - mol2.SetProp("Converged", "true") - mol2.SetProp("E_tot", "-10.5") # |dE| >> tol - - unique_mols = filter_unique([mol, mol2], crit=0.3) - assert len(unique_mols) == 2 - - def test_missing_energy_falls_back_to_rmsd_only(self): - """Without E_tot the energy guard cannot apply -> RMSD-only dedup. - - Preserves the legacy behavior for callers that do not set E_tot. - """ - - mol = Chem.MolFromSmiles("CCO") - mol = Chem.AddHs(mol) - AllChem.EmbedMolecule(mol, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol) - mol.SetProp("Converged", "true") # no E_tot set - - mol2 = Chem.Mol(mol) - mol2.SetProp("Converged", "true") - - unique_mols = filter_unique([mol, mol2], crit=0.3) - assert len(unique_mols) == 1 - - def test_filter_different_conformers(self): - """Test that different conformers are kept.""" - - mol1 = Chem.MolFromSmiles("CCCCCC") # Hexane - flexible - mol1 = Chem.AddHs(mol1) - AllChem.EmbedMolecule(mol1, randomSeed=42) - mol1.SetProp("Converged", "true") - - mol2 = Chem.MolFromSmiles("CCCCCC") - mol2 = Chem.AddHs(mol2) - AllChem.EmbedMolecule(mol2, randomSeed=123) - mol2.SetProp("Converged", "true") - - # Generate very different conformers by using different seeds - # and moving atoms around - conf = mol2.GetConformer() - pos = conf.GetAtomPosition(0) - conf.SetAtomPosition(0, (pos.x + 0.5, pos.y, pos.z)) - - mols = [mol1, mol2] - unique_mols = filter_unique(mols, crit=0.3) - - # Should keep both (or at least not crash) - assert len(unique_mols) >= 1 - - def test_two_diastereomers_are_never_merged(self): - """A distinct compound must survive dedup, however close its geometry. - - The same guarantee ``tests/test_filtering.py`` asserts for - ``filter_unique_optimized``, asserted here because this is the other - duplicate filter and it applies the identical RMSD-plus-energy criterion. - Fixing one path and not the other would leave the defect reachable - through ``ConformerRanker(use_optimized_filtering=False)``. - - cis/trans-4-tert-butylcyclohexanol: heavy-atom RMSD between the two - diastereomers was measured at 0.300 A, i.e. at the 0.3 A default - threshold. ``crit`` is opened wide here so RMSD and energy both say - "duplicate" and the only thing that can keep the pair apart is the fact - that they are different compounds. - """ - from Auto3D.utils.energy import set_e_tot_from_ev - - def build(smiles: str) -> Chem.Mol: - mol = Chem.AddHs(Chem.MolFromSmiles(smiles)) - AllChem.EmbedMolecule(mol, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol) - mol.SetProp("Converged", "true") - set_e_tot_from_ev(mol, -10.0) # identical energies - return mol - - cis = build("O[C@H]1CC[C@@H](CC1)C(C)(C)C") - trans = build("O[C@H]1CC[C@H](CC1)C(C)(C)C") - assert Chem.MolToSmiles(cis) != Chem.MolToSmiles(trans), "test premise" - - assert len(filter_unique([cis, trans], crit=10.0)) == 2, ( - "the legacy filter merged two distinct diastereomers, so an input " - "molecule vanished from the output with no record" - ) - - def test_duplicate_conformers_of_one_stereoisomer_still_collapse(self): - """The other half: the stereo guard must narrow dedup, not disable it.""" - from Auto3D.utils.energy import set_e_tot_from_ev - - def build() -> Chem.Mol: - mol = Chem.AddHs(Chem.MolFromSmiles("O[C@H]1CC[C@@H](CC1)C(C)(C)C")) - AllChem.EmbedMolecule(mol, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol) - mol.SetProp("Converged", "true") - set_e_tot_from_ev(mol, -10.0) - return mol - - assert len(filter_unique([build(), build()], crit=10.0)) == 1, ( - "duplicate conformers of one stereoisomer survived, so the stereo " - "guard has switched dedup off rather than narrowing it" - ) - - def test_filter_unconverged_removed(self): - """Test that unconverged structures are removed.""" - - mol1 = Chem.MolFromSmiles("CCO") - mol1 = Chem.AddHs(mol1) - AllChem.EmbedMolecule(mol1, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol1) - mol1.SetProp("Converged", "true") - - mol2 = Chem.MolFromSmiles("CCO") - mol2 = Chem.AddHs(mol2) - AllChem.EmbedMolecule(mol2, randomSeed=123) - mol2.SetProp("Converged", "false") # Not converged - - mols = [mol1, mol2] - unique_mols = filter_unique(mols, crit=0.3) - - # Only converged one should remain - assert len(unique_mols) == 1 - assert unique_mols[0].GetProp("Converged").lower() == "true" - - def test_filter_empty_list(self): - """Test filtering empty list returns empty list.""" - - unique_mols = filter_unique([], crit=0.3) - assert len(unique_mols) == 0 - - def test_filter_custom_threshold(self): - """Test that custom RMSD threshold works.""" - - mol1 = Chem.MolFromSmiles("CCO") - mol1 = Chem.AddHs(mol1) - AllChem.EmbedMolecule(mol1, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol1) - mol1.SetProp("Converged", "true") - - mol2 = Chem.Mol(mol1) - mol2.SetProp("Converged", "true") - - mols = [mol1, mol2] - - # With very small threshold, might keep both - unique_mols_small = filter_unique(mols, crit=0.0001) - # With large threshold, definitely keep only one - unique_mols_large = filter_unique(mols, crit=10.0) - - # Large threshold should definitely merge identical mols - assert len(unique_mols_large) == 1 - # A tighter threshold can never merge MORE than a looser one -- the - # discarded half of this test's own computation, now actually checked. - assert len(unique_mols_small) >= len(unique_mols_large) - - def test_filter_unique_removehs_is_linear_and_nondestructive(self, monkeypatch): - """Legacy filter_unique strips Hs once per molecule (not per comparison) and - returns the originals with explicit H + exact positions intact.""" - import numpy as np - from rdkit import Chem - from rdkit.Chem import AllChem - - from Auto3D import filtering - - base = Chem.AddHs(Chem.MolFromSmiles("CCCCO")) - cids = AllChem.EmbedMultipleConfs(base, numConfs=5, randomSeed=1) - mols = [] - for cid in cids: - m = Chem.Mol(base, confId=int(cid)) - m.SetProp("Converged", "true") - mols.append(m) - n_atoms = base.GetNumAtoms() - orig_pos = {id(m): m.GetConformer().GetPositions().copy() for m in mols} - - calls = {"n": 0} - real_removehs = filtering.Chem.RemoveHs - - def counting(mol, *a, **k): - calls["n"] += 1 - return real_removehs(mol, *a, **k) - - monkeypatch.setattr(filtering.Chem, "RemoveHs", counting) - - result = filtering.filter_unique(mols, crit=0.01) - assert calls["n"] == len(mols) # once per input, never per pair - assert len(result) == len(mols) - for m in result: - assert m.GetNumAtoms() == n_atoms - assert any(a.GetAtomicNum() == 1 for a in m.GetAtoms()) - assert np.array_equal(m.GetConformer().GetPositions(), orig_pos[id(m)]) - - def test_rmsd_failure_keeps_both(self, monkeypatch): - """An incomparable pair (RMSD raises) must NOT be treated as a duplicate. - - When GetBestRMS raises RuntimeError, filter_unique must treat the pair - as distinct (rmsd = inf) and keep both, mirroring the fix already in - filtering._filter_within_cluster. The previous behavior (rmsd = 0) - made distinct conformers look like perfect duplicates and dropped one. - """ - from Auto3D import filtering - - def make(name): - m = Chem.AddHs(Chem.MolFromSmiles("CCO")) - AllChem.EmbedMolecule(m, randomSeed=abs(hash(name)) % 1000) - AllChem.MMFFOptimizeMolecule(m) - m.SetProp("_Name", name) - m.SetProp("Converged", "true") - return m - - def boom(*args, **kwargs): - raise RuntimeError("GetBestRMS failed") - - # filter_unique calls rdMolAlign.GetBestRMS via the filtering module. - monkeypatch.setattr(filtering.rdMolAlign, "GetBestRMS", boom) - - mols = [make("a"), make("b")] - unique_mols = filtering.filter_unique(mols, crit=0.3) - assert len(unique_mols) == 2 # incomparable pair must NOT be dropped diff --git a/tests/test_filtering.py b/tests/test_filtering.py index 6974ca7c..759858d0 100644 --- a/tests/test_filtering.py +++ b/tests/test_filtering.py @@ -1,12 +1,20 @@ #!/usr/bin/env python -"""Tests for optimized RMSD filtering with energy clustering.""" +"""Tests for the single conformer filter (RMSD dedup with energy clustering).""" from __future__ import annotations +import os + import pytest from rdkit import Chem from rdkit.Chem import AllChem -from Auto3D.filtering import filter_unique_optimized, _filter_within_cluster +from Auto3D.filtering import ( + DROP_REASONS, + FilterResult, + _filter_within_cluster, + filter_conformers, + filter_unique_optimized, +) from Auto3D.utils.energy import set_e_tot_from_ev @@ -331,38 +339,43 @@ def test_small_energy_window_creates_separate_clusters(self): assert len(result) == 2 +def _energyless(smiles: str, seed: int = 42, name: str = "") -> Chem.Mol: + """An embedded, converged conformer carrying NO 'E_tot' property. + + This is what an SDF Auto3D's optimizer did not write can look like: a + hand-built conformer set, or an export that names its energy field + something else. ``ConformerRanker`` refuses such a record up front + (``InputValidationError``), so this shape only reaches the filter through + a direct API call -- which is exactly the caller the two filters used to + disagree for. + """ + mol = Chem.AddHs(Chem.MolFromSmiles(smiles)) + AllChem.EmbedMolecule(mol, randomSeed=seed) + AllChem.MMFFOptimizeMolecule(mol) + mol.SetProp("Converged", "true") + if name: + mol.SetProp("_Name", name) + assert not mol.HasProp("E_tot"), "helper premise: no energy property" + return mol + + class TestMissingEnergyPropertyMustNotCrash: - """filter_unique_optimized must tolerate a record with no 'E_tot', the - way the legacy ``filtering.filter_unique`` already does. - - KNOWN DEFECT (found during cluster E brainstorming, not fixed by this - lane): ``filtering.py:75`` sorts the valid-mols list by - ``Auto3D.utils.energy.e_tot_ev``, which RAISES (KeyError/ValueError) for - a molecule with no usable 'E_tot' property. ``_filter_within_cluster``'s - own energy guard, two dozen lines later in the same module, instead uses - the tolerant ``try_e_tot_ev`` and treats a missing energy as "fall back - to RMSD only". ``filtering.filter_unique`` (the OTHER conformer - filter, sharing the same duplicate criterion since 4.0.1) also uses - ``try_e_tot_ev`` throughout and does not crash on this input. So the two - filters diverge on malformed input: the same list of mols that - ``filter_unique`` happily filters crashes ``filter_unique_optimized``. - - This matters here specifically because cluster B5 is about to delete one - of the two filters, and the survivor is the stricter (crashing) one -- - fixing filtering.py's sort key to use ``try_e_tot_ev``, matching its own - energy guard and the legacy filter, is what should make this pass. + """The one conformer filter must tolerate a record with no 'E_tot'. + + ``filtering.py`` used to sort the valid-mols list by + ``Auto3D.utils.energy.e_tot_ev``, which RAISES (KeyError/ValueError) for a + molecule with no usable 'E_tot'. ``_filter_within_cluster``'s own energy + guard, two dozen lines later in the same module, instead used the tolerant + ``try_e_tot_ev`` and treated a missing energy as "fall back to RMSD only", + as did the legacy all-pairs ``filter_unique`` throughout. So the same list + of mols one filter happily filtered crashed the other -- and the survivor + of the two was the crashing one. + + Tolerating the record must not mean *inventing* an energy for it: a + missing 'E_tot' read as 0.0 would sort a garbage record to the front of a + list of negative energies and hand it to ``top_k`` as the global minimum. """ - @pytest.mark.xfail( - strict=True, - reason=( - "filtering.py:75 sorts by e_tot_ev (raises KeyError for a mol " - "with no 'E_tot' property) instead of the tolerant try_e_tot_ev " - "that _filter_within_cluster's own energy guard and the legacy " - "filtering.filter_unique both use -- the two conformer " - "filters disagree on malformed input (cluster E brainstorm defect)." - ), - ) def test_missing_e_tot_property_does_not_crash(self): mol_no_energy = Chem.AddHs(Chem.MolFromSmiles("CCO")) AllChem.EmbedMolecule(mol_no_energy, randomSeed=42) @@ -380,20 +393,116 @@ def test_missing_e_tot_property_does_not_crash(self): # cannot be deduped by energy -- it must survive alongside the other. assert len(result) == 2 + def test_a_garbled_e_tot_is_tolerated_the_same_way(self): + """``try_e_tot_ev`` swallows ValueError too, so a non-numeric property + must take the same path as an absent one rather than crash.""" + mol_garbled = _energyless("CCO", seed=42) + mol_garbled.SetProp("E_tot", "not-a-number") + + mol_with_energy = _create_mol_with_energy("CC", -10.0) + + result = filter_unique_optimized( + [mol_garbled, mol_with_energy], rmsd_threshold=0.3 + ) + assert len(result) == 2 + + def test_an_energyless_record_sorts_last_whatever_the_real_energies_are(self): + """A tolerated missing energy must not be *compared* as an energy. + + The tolerant path needs some placeholder to sort on; the danger is that + the placeholder takes part in the comparison. A record with no ``E_tot`` + that sorts as if its energy were 0.0 lands ahead of every genuine + structure whose energy is above 0.0 -- and the first element of the + filter's output is what ``top_k``/``top_window`` treat as the global + minimum: the reference ``E_rel`` is measured from, and the single + structure a ``k=1`` request returns. + + ``E_tot`` is a user-supplied SD property, so "every real energy is a + large negative number" is a convention, not a guarantee -- hence the + deliberately positive energy below. The record with no energy must sort + last regardless of where the real energies fall relative to the + placeholder. + """ + low = _create_mol_with_energy("CCO", -100.0) + high = _create_mol_with_energy("CCCO", -50.0) + positive = _create_mol_with_energy("CCCCCCO", +25.0) + unknown = _energyless("CCCCO", seed=7) + + result = filter_unique_optimized( + [unknown, positive, high, low], rmsd_threshold=0.3 + ) + + assert len(result) == 4 + assert [Chem.MolToSmiles(Chem.RemoveHs(m)) for m in result] == [ + Chem.MolToSmiles(Chem.RemoveHs(low)), + Chem.MolToSmiles(Chem.RemoveHs(high)), + Chem.MolToSmiles(Chem.RemoveHs(positive)), + Chem.MolToSmiles(Chem.RemoveHs(unknown)), + ], "the record with no energy must sort after every record that has one" + + def test_two_energyless_duplicates_still_collapse(self): + """The inverse assertion, and the reason the two above are safe. + + A filter that declared every energy-less record distinct would pass + both tests above while silently switching duplicate removal off for + this whole class of input. Two identical geometries with no energy at + all must still collapse to one: the energy guard cannot apply, so RMSD + alone decides -- which is precisely what the legacy all-pairs filter + did. + """ + first = _energyless("CCO", seed=42) + second = _energyless("CCO", seed=42) + + result = filter_unique_optimized([first, second], rmsd_threshold=0.3) -class TestFilterUniqueBehavior: - """Tests verifying behavior matches original filter_unique.""" + assert len(result) == 1, ( + "two bit-identical energy-less conformers both survived, so " + "tolerating a missing energy has disabled dedup for this input " + "rather than falling back to RMSD only" + ) - def test_matches_original_for_simple_case(self): - """Should produce same results as original filter_unique for basic cases. + def test_an_energyless_record_is_compared_against_one_that_has_energy(self): + """Energy-less records are compared against EVERYTHING, not just each + other. - Both filters operate on conformers of the *same* molecule (that is the - real Auto3D contract). Use genuinely distinct conformers of one molecule - so RMSD is well-defined and comparable across both implementations. + No energy gap can prove a pair is not a duplicate when one side has no + energy, so the partitioning that makes the filter sub-quadratic has no + licence to separate such a record from anything. The legacy all-pairs + filter compared it to every survivor; so must this one. """ - from Auto3D.filtering import filter_unique + with_energy = _create_mol_with_energy("CCO", -100.0) + without = _energyless("CCO", seed=42) + # Same compound, same embedding seed and force field -> RMSD ~= 0. + assert Chem.MolToSmiles(Chem.RemoveHs(with_energy)) == Chem.MolToSmiles( + Chem.RemoveHs(without) + ), "test premise: same compound" - def conformer(seed: float, energy_ev: float) -> Chem.Mol: + result = filter_unique_optimized( + [with_energy, without], rmsd_threshold=0.3, energy_cluster_window=0.01 + ) + + assert len(result) == 1, ( + "a record with no energy escaped comparison against a duplicate " + "that has one, because the energy partitioning separated them" + ) + + +class TestTheSurvivingFilterKeepsTheLegacyVerdicts: + """Values recorded from the legacy all-pairs ``filter_unique`` before it was + deleted (cluster B5 phase 4a). + + Auto3D carried two conformer filters with the same duplicate criterion, each + acting as the other's oracle, until 4.1.0. These cases were run against BOTH + and are asserted here as literals so the surviving filter's verdicts stay + pinned now that there is nothing left to compare against. + """ + + def test_distinct_conformers_of_one_molecule_all_survive(self): + """Three genuinely different conformers of a flexible chain: 3 kept. + + Both filters returned 3 for this input. + """ + def conformer(seed: int, energy_ev: float) -> Chem.Mol: m = Chem.AddHs(Chem.MolFromSmiles("CCCCCCO")) # flexible chain AllChem.EmbedMolecule(m, randomSeed=seed) AllChem.MMFFOptimizeMolecule(m) @@ -403,15 +512,121 @@ def conformer(seed: float, energy_ev: float) -> Chem.Mol: mols = [conformer(42, -12.0), conformer(7, -11.0), conformer(123, -10.0)] - original_result = filter_unique(mols, crit=0.3) - optimized_result = filter_unique_optimized( - mols, - rmsd_threshold=0.3, - energy_cluster_window=100.0 # Large window = single cluster = same behavior + assert len(filter_unique_optimized(mols, rmsd_threshold=0.3)) == 3 + # A single cluster (huge window) must give the same answer -- that + # equivalence is the whole justification for the energy partitioning. + assert len( + filter_unique_optimized( + mols, rmsd_threshold=0.3, energy_cluster_window=100.0 + ) + ) == 3 + + def test_a_malformed_mixed_list_keeps_the_recorded_three(self): + """The input the two filters used to DISAGREE about. + + An energy-bearing conformer, an energy-less duplicate of it, an + energy-less record of a different compound, and a distinct third + compound. The legacy filter kept three of the four -- merging the + energy-less ethanol into the one that has an energy -- and so must this + one. + """ + a = _create_mol_with_energy("CCO", -100.0) + a.SetProp("_Name", "ethanol_with_energy") + b = _energyless("CCO", seed=42, name="ethanol_no_energy") + c = _energyless("CCCCO", seed=11, name="butanol_no_energy") + d = _create_mol_with_energy("CCCCCCO", -80.0) + d.SetProp("_Name", "heptanol_with_energy") + + kept = filter_unique_optimized([a, b, c, d], rmsd_threshold=0.3) + + assert {m.GetProp("_Name") for m in kept} == { + "ethanol_with_energy", + "butanol_no_energy", + "heptanol_with_energy", + } + + def test_an_rmsd_threshold_sweep_straddling_a_measured_pair(self): + """A threshold below the pair's actual RMSD keeps both; above, merges. + + Recorded from both filters. Sweeping across the *measured* RMSD -- not + a fixed number -- is what makes this fail for a filter that ignores + ``rmsd_threshold`` altogether, which an equal-length comparison at a + single threshold would not. + """ + from rdkit.Chem import rdMolAlign + + first = _create_mol_with_energy("CCCCCCCC", -100.0) # octane + second = Chem.Mol(first) + AllChem.EmbedMolecule(second, randomSeed=99) + AllChem.MMFFOptimizeMolecule(second) + set_e_tot_from_ev(second, -100.0) # energies agree -> RMSD decides + second.SetProp("Converged", "true") + + rmsd = rdMolAlign.GetBestRMS(Chem.RemoveHs(first), Chem.RemoveHs(second)) + assert rmsd > 0.05, "test premise: the conformers must be distinct" + + for crit, expected in ((rmsd / 2, 2), (rmsd * 2, 1)): + kept = filter_unique_optimized([first, second], rmsd_threshold=crit) + assert len(kept) == expected, f"at rmsd_threshold={crit}" + + def test_two_energyless_conformers_dedup_by_rmsd_alone(self): + """Ported from the legacy filter's own suite. + + Neither record carries ``E_tot``, so the energy half of the duplicate + criterion cannot apply and RMSD alone decides. The legacy filter kept + one; so must this one. + """ + mol = _energyless("CCO", seed=42) + duplicate = Chem.Mol(mol) + duplicate.SetProp("Converged", "true") + + assert len(filter_unique_optimized([mol, duplicate], rmsd_threshold=0.3)) == 1 + + +class TestConvergencePropertyAbsenceFiltersLikeTrue: + """A whole SDF that carries no 'Converged' property must filter exactly as + the same records marked Converged=True do. + + Ported from the legacy filter's suite (it lived beside the validation tests + because ``filter_unique`` used to live in ``utils/validation.py``). Only + ``batchopt`` writes ``Converged``; an ``opt_geometry`` output, an + ORCA/Gaussian export or a hand-built conformer set carries none, and + treating that as "did not converge" deleted every record. + """ + + _SDF = os.path.join( + os.path.dirname(os.path.abspath(__file__)), "files", "example.sdf" + ) + + def _mols(self) -> list[Chem.Mol]: + supp = Chem.SDMolSupplier(self._SDF, removeHs=False) + return [mol for mol in supp if mol is not None] + + def test_an_unflagged_file_keeps_what_a_flagged_one_keeps(self): + flagged = self._mols() + for mol in flagged: + mol.SetProp("Converged", "True") + expected = len(filter_unique_optimized(flagged, rmsd_threshold=0.3)) + assert expected >= 1, "test premise: the flagged file must keep something" + + unflagged = self._mols() + for mol in unflagged: + mol.ClearProp("Converged") + assert not mol.HasProp("Converged") + + result = filter_unique_optimized(unflagged, rmsd_threshold=0.3) + assert len(result) == expected, ( + f"{len(unflagged)} record(s) with no 'Converged' property kept " + f"{len(result)}, but the same records marked Converged=True keep " + f"{expected}" ) - # Same molecule, well-defined RMSD: both implementations must agree. - assert len(original_result) == len(optimized_result) + def test_an_explicit_false_still_empties_the_selection(self): + """The inverse: absence is not failure, but a stated failure is.""" + mols = self._mols() + for mol in mols: + mol.SetProp("Converged", "False") + assert filter_unique_optimized(mols, rmsd_threshold=0.3) == [] def test_filter_within_cluster_removehs_is_linear_and_nondestructive(monkeypatch): @@ -482,3 +697,101 @@ def boom(*a, **k): cluster = [make("a", -1.0), make("b", -0.9)] kept = _filter_within_cluster(cluster, rmsd_threshold=0.3) assert len(kept) == 2 # incomparable pair must NOT be dropped + + +class TestTheFilterSaysWhyItDroppedThings: + """``filter_conformers`` reports a count per reason, not just a survivor list. + + Returning a bare list is what let ``ranking`` tell a user "No structure + converged" for a species whose conformers were every one of them dropped + for *stereochemistry* -- a message that points at ``--opt-steps`` and + ``--convergence-threshold`` for a problem neither can fix. + """ + + def test_each_reason_is_counted_under_its_own_name(self): + good = _create_mol_with_energy("CCO", -100.0) + unconverged = _create_mol_with_energy("CCO", -99.0, converged=False) + stereo_changed = _create_mol_with_energy("CCO", -98.0) + stereo_changed.SetProp("Stereo_changed", "true") + broken = _create_mol_with_energy("CC", -97.0) + conf = broken.GetConformer() + pos = conf.GetAtomPosition(0) + conf.SetAtomPosition(0, (pos.x + 5.0, pos.y, pos.z)) + duplicate = _create_mol_with_energy("CCO", -100.0) # same as `good` + + result = filter_conformers( + [None, good, unconverged, stereo_changed, broken, duplicate], + rmsd_threshold=0.3, + ) + + assert [m is good for m in result.kept] == [True] + assert result.dropped == { + "unparsed": 1, + "unconverged": 1, + "stereochemistry": 1, + "connectivity": 1, + "duplicate": 1, + } + + def test_nothing_dropped_reports_nothing(self): + """The inverse: a clean input must not manufacture a reason. + + A result object that always carried a non-empty ``dropped`` would make + every ranking message name a cause that did not happen. + """ + result = filter_conformers( + [_create_mol_with_energy("CCO", -100.0)], rmsd_threshold=0.3 + ) + assert result.dropped == {} + assert result.reasons == () + assert result.summary() == "" + + def test_summary_names_every_reason_that_fired_in_declared_order(self): + result = FilterResult( + kept=[], + dropped={"duplicate": 3, "unconverged": 2, "stereochemistry": 1}, + ) + assert result.reasons == ("unconverged", "stereochemistry", "duplicate") + assert result.summary() == ( + "2 marked Converged=false, " + "1 changed stereochemistry during optimization, " + "3 duplicates of a kept conformer" + ) + + def test_zero_counts_are_not_reported_as_reasons(self): + result = FilterResult(kept=[], dropped={"unconverged": 0, "duplicate": 2}) + assert result.reasons == ("duplicate",) + assert result.summary() == "2 duplicates of a kept conformer" + + def test_an_unknown_reason_is_refused_at_construction(self): + """DROP_REASONS is the vocabulary, enforced. + + Without this, a producer misspelling a reason contributes a drop that + ``summary()`` silently omits, so a user is told fewer conformers went + missing than actually did. + """ + with pytest.raises(ValueError, match="unknown filter drop reason"): + FilterResult(kept=[], dropped={"unconvrged": 1}) + + def test_every_declared_reason_has_a_phrase(self): + """A reason with no phrase would raise KeyError from inside summary() + -- while reporting a diagnostic, which is the worst place to fail.""" + for reason in DROP_REASONS: + assert FilterResult(kept=[], dropped={reason: 1}).summary() + + def test_truncation_is_not_a_drop_reason(self): + """`k` cutting the list short is selection, not a filter drop. + + Nothing is missing to explain: those conformers are valid and unique, + they lost the ranking. Counting them would make every ``k=1`` run report + drops it should not. + """ + assert "truncated" not in DROP_REASONS + distinct = [ + _create_mol_with_energy("CCO", -100.0), + _create_mol_with_energy("CCCO", -90.0), + _create_mol_with_energy("CCCCO", -80.0), + ] + result = filter_conformers(distinct, rmsd_threshold=0.3) + assert len(result.kept) == 3 + assert result.dropped == {} diff --git a/tests/test_id_mapping.py b/tests/test_id_mapping.py index 313fe297..800086c2 100644 --- a/tests/test_id_mapping.py +++ b/tests/test_id_mapping.py @@ -133,6 +133,14 @@ def test_encode_ids_refuses_to_overwrite_an_existing_file(self, tmp_path): `out_dir` below), but this check keeps the guarantee attached to the function itself, so a caller taking the default location cannot reintroduce the defect. + + The refusal now comes from the shared + `Auto3D.utils.output_guard.check_output_overwrite` (hence the + "already exists" wording) rather than a bespoke message here, so + `encode_ids`, `decode_ids` and `tautomer.select_tautomers` state the + same policy the same way -- and `overwrite=True` can lift it, which the + unconditional refusal this replaced offered no way to do. + `tests/test_output_overwrite_gates.py` covers both directions. """ from Auto3D.exceptions import ConfigurationError @@ -141,7 +149,7 @@ def test_encode_ids_refuses_to_overwrite_an_existing_file(self, tmp_path): users_file = tmp_path / "mols_encoded.smi" users_file.write_bytes(b"IRREPLACEABLE USER DATA\n") - with pytest.raises(ConfigurationError, match="would overwrite"): + with pytest.raises(ConfigurationError, match="already exists"): encode_ids(str(p)) assert users_file.read_bytes() == b"IRREPLACEABLE USER DATA\n" diff --git a/tests/test_output_overwrite_gates.py b/tests/test_output_overwrite_gates.py new file mode 100644 index 00000000..f6ee127c --- /dev/null +++ b/tests/test_output_overwrite_gates.py @@ -0,0 +1,259 @@ +#!/usr/bin/env python +"""Overwrite gates on the public writers that derive their own output names. + +``Chem.SDWriter(path)`` and ``open(path, "w")`` truncate on open, and four +public functions used to do that with no consent gate at all: + +* ``Auto3D.tautomer.select_tautomers`` -- ``select_tautomers("/data/results.sdf", + k=1)`` replaced ``/data/results_top_tautomers.sdf``, a name it invented, with + this call's selection. +* ``Auto3D.id_mapping.decode_ids`` -- same shape, for ``_out.sdf``. +* ``Auto3D.utils.smi_io.smiles2smi`` -- the caller named the file here, so the + gate defaults *open*; what matters is that it can be closed. +* ``Auto3D.id_mapping.encode_ids`` -- refused unconditionally, with no way for a + caller to say yes. + +The policy: permissive where the caller named the file, restrictive where Auto3D +invented the name. Every gate is keyword-only, and every one is asserted in both +directions -- a gate that refused legitimate writes would satisfy the "refuses" +half of each pair and break the pipeline. +""" +from __future__ import annotations + +import pytest +from rdkit import Chem +from rdkit.Chem import AllChem + +from Auto3D.exceptions import ConfigurationError +from Auto3D.id_mapping import decode_ids, encode_ids +from Auto3D.tautomer import select_tautomers +from Auto3D.utils.energy import set_e_tot_from_ev +from Auto3D.utils.smi_io import smiles2smi + +PRECIOUS = b"IRREPLACEABLE USER DATA\n" + + +def _tautomer_sdf(path, names=("mol@taut1", "mol@taut2")) -> str: + """An SDF shaped the way ``select_tautomers`` expects: ``id@tautN`` names.""" + with Chem.SDWriter(str(path)) as writer: + for i, name in enumerate(names): + mol = Chem.AddHs(Chem.MolFromSmiles("CCO")) + AllChem.EmbedMolecule(mol, randomSeed=42 + i) + mol.SetProp("_Name", name) + set_e_tot_from_ev(mol, -10.0 - i) + writer.write(mol) + return str(path) + + +def _encoded_sdf(path) -> str: + """An SDF shaped the way ``decode_ids`` expects: numeric name + ID prop.""" + with Chem.SDWriter(str(path)) as writer: + mol = Chem.AddHs(Chem.MolFromSmiles("CCO")) + AllChem.EmbedMolecule(mol, randomSeed=42) + mol.SetProp("_Name", "0") + mol.SetProp("ID", "0_0_0") + writer.write(mol) + return str(path) + + +class TestSelectTautomers: + """Auto3D invents ``_top_tautomers.sdf``, so the gate defaults shut.""" + + def test_it_refuses_to_replace_an_existing_top_tautomers_file(self, tmp_path): + """The concrete hazard from the follow-ups list.""" + sdf = _tautomer_sdf(tmp_path / "results.sdf") + derived = tmp_path / "results_top_tautomers.sdf" + derived.write_bytes(PRECIOUS) + + with pytest.raises(ConfigurationError, match="already exists"): + select_tautomers(sdf, k=1) + + assert derived.read_bytes() == PRECIOUS + + def test_it_writes_when_the_derived_name_is_free(self, tmp_path): + """Inverse: the ordinary case must still go through untouched. + + Every route Auto3D itself takes lands here -- ``get_stable_tautomers`` + passes ``main()``'s output from a job directory created fresh for that + run -- so a gate that refused unconditionally would break the whole + tautomer pipeline while passing the test above. + """ + sdf = _tautomer_sdf(tmp_path / "results.sdf") + out = select_tautomers(sdf, k=1) + + assert out == str(tmp_path / "results_top_tautomers.sdf") + assert len( + [m for m in Chem.SDMolSupplier(out, removeHs=False) if m is not None] + ) == 1 + + def test_overwrite_true_replaces_it(self, tmp_path): + """Inverse: the gate is a consent gate the caller can lift.""" + sdf = _tautomer_sdf(tmp_path / "results.sdf") + derived = tmp_path / "results_top_tautomers.sdf" + derived.write_bytes(PRECIOUS) + + out = select_tautomers(sdf, k=1, overwrite=True) + + assert out == str(derived) + assert derived.read_bytes() != PRECIOUS + + def test_overwrite_is_keyword_only(self, tmp_path): + """``select_tautomers(sdf, k, window)`` is the documented positional + signature; a fourth positional must not silently become ``overwrite``.""" + sdf = _tautomer_sdf(tmp_path / "results.sdf") + with pytest.raises(TypeError): + select_tautomers(sdf, 1, None, True) + + +class TestDecodeIds: + """Auto3D invents ``_out.sdf``, so the gate defaults shut.""" + + def test_it_refuses_to_replace_an_existing_out_file(self, tmp_path): + sdf = _encoded_sdf(tmp_path / "mols_3d_encoded.sdf") + derived = tmp_path / "mols_out.sdf" + derived.write_bytes(PRECIOUS) + + with pytest.raises(ConfigurationError, match="already exists"): + decode_ids(sdf, {"mol_a": 0}) + + assert derived.read_bytes() == PRECIOUS + + def test_it_writes_when_the_derived_name_is_free(self, tmp_path): + """Inverse: this is the call ``WorkflowOrchestrator`` makes on every run. + + It writes into the job directory it created with a bare ``mkdir()`` + moments earlier, so the derived name is always free there -- but only if + the gate permits a free name. + """ + sdf = _encoded_sdf(tmp_path / "mols_3d_encoded.sdf") + out = decode_ids(sdf, {"mol_a": 0}) + + assert out == str(tmp_path / "mols_out.sdf") + written = [ + m for m in Chem.SDMolSupplier(out, removeHs=False) if m is not None + ] + assert [m.GetProp("_Name") for m in written] == ["mol_a"] + + def test_overwrite_true_replaces_it(self, tmp_path): + sdf = _encoded_sdf(tmp_path / "mols_3d_encoded.sdf") + derived = tmp_path / "mols_out.sdf" + derived.write_bytes(PRECIOUS) + + out = decode_ids(sdf, {"mol_a": 0}, overwrite=True) + + assert out == str(derived) + assert derived.read_bytes() != PRECIOUS + + def test_overwrite_is_keyword_only(self, tmp_path): + sdf = _encoded_sdf(tmp_path / "mols_3d_encoded.sdf") + with pytest.raises(TypeError): + decode_ids(sdf, {"mol_a": 0}, True) + + +class TestEncodeIds: + """Already refused unconditionally; now it refuses *consistently*. + + Same keyword, same default, same exception and same message as the other + three -- and, unlike before, a caller who means it can say so. + """ + + def test_it_still_refuses_by_default(self, tmp_path): + smi = tmp_path / "mols.smi" + smi.write_text("CCO a\n") + derived = tmp_path / "mols_encoded.smi" + derived.write_bytes(PRECIOUS) + + with pytest.raises(ConfigurationError, match="already exists"): + encode_ids(str(smi)) + + assert derived.read_bytes() == PRECIOUS + + def test_overwrite_true_replaces_it(self, tmp_path): + """The half that was impossible before: there was no way to consent.""" + smi = tmp_path / "mols.smi" + smi.write_text("CCO a\n") + derived = tmp_path / "mols_encoded.smi" + derived.write_bytes(PRECIOUS) + + new_path, mapping = encode_ids(str(smi), overwrite=True) + + assert new_path == str(derived) + assert mapping == {"a": 0} + assert derived.read_text() == "CCO 0\n" + + def test_it_writes_when_the_derived_name_is_free(self, tmp_path): + """Inverse: the pipeline's own call, which must keep working.""" + smi = tmp_path / "mols.smi" + smi.write_text("CCO a\nCCC b\n") + new_path, mapping = encode_ids(str(smi)) + assert mapping == {"a": 0, "b": 1} + assert new_path == str(tmp_path / "mols_encoded.smi") + + def test_overwrite_is_keyword_only(self, tmp_path): + smi = tmp_path / "mols.smi" + smi.write_text("CCO a\n") + with pytest.raises(TypeError): + encode_ids(str(smi), None, True) + + +class TestSmiles2Smi: + """The caller named this file, so the gate defaults OPEN. + + ``smiles2mols`` writes it into a ``TemporaryDirectory`` on every call, and a + caller who passes an explicit path has already chosen it. Defaulting shut + here would be a gate on a decision the caller already made. + """ + + def test_the_default_still_overwrites(self, tmp_path): + """Inverse first, because this is the back-compat guarantee.""" + out = tmp_path / "mols.smi" + out.write_bytes(PRECIOUS) + + assert smiles2smi(["CCO"], str(out)) == str(out) + assert out.read_bytes() != PRECIOUS + assert out.read_text().startswith("CCO ") + + def test_overwrite_false_refuses(self, tmp_path): + out = tmp_path / "mols.smi" + out.write_bytes(PRECIOUS) + + with pytest.raises(ConfigurationError, match="already exists"): + smiles2smi(["CCO"], str(out), overwrite=False) + + assert out.read_bytes() == PRECIOUS + + def test_overwrite_false_still_writes_a_free_path(self, tmp_path): + out = tmp_path / "mols.smi" + assert smiles2smi(["CCO"], str(out), overwrite=False) == str(out) + assert out.read_text().startswith("CCO ") + + def test_overwrite_is_keyword_only(self, tmp_path): + with pytest.raises(TypeError): + smiles2smi(["CCO"], str(tmp_path / "mols.smi"), False) + + +class TestTheGateRefusesBeforeDoingTheWork: + """A gate applied just before the write is a gate that wasted the run. + + ``select_tautomers`` reads and groups the whole input before it opens the + writer; refusing only there means the user waits for the work, then loses + it. The check must come first. + """ + + def test_select_tautomers_checks_before_reading_the_input( + self, tmp_path, monkeypatch + ): + import Auto3D.tautomer as tautomer + + sdf = _tautomer_sdf(tmp_path / "results.sdf") + (tmp_path / "results_top_tautomers.sdf").write_bytes(PRECIOUS) + + def _never(*args, **kwargs): + raise AssertionError( + "select_tautomers read its input before checking the output path" + ) + + monkeypatch.setattr(tautomer.Chem, "SDMolSupplier", _never) + + with pytest.raises(ConfigurationError, match="already exists"): + select_tautomers(sdf, k=1) diff --git a/tests/test_padding_invariance.py b/tests/test_padding_invariance.py index a96eb50b..cbe54d4a 100644 --- a/tests/test_padding_invariance.py +++ b/tests/test_padding_invariance.py @@ -2,9 +2,11 @@ Every engine masks padded slots differently: AIMNet2 pads species with 0 and relies on Z=0 being unused, while ANI2x/ANI2xt pad with -1 and rely on that -being torchani's masked-atom sentinel. batch_opt/ANI2xt_no_rep.py:167-172 -documents that second assumption as unverified. These tests verify it (audit -M32, C13). +being torchani's masked-atom sentinel. ``ANI2xt.forward`` in +``batch_opt/ANI2xt_no_rep.py`` documents that second assumption (species +== -1 surviving the periodic-table remap unchanged and relying on +TorchANI's masked-atom convention) as depended-upon but not independently +verified there. These tests verify it (audit M32, C13). """ from __future__ import annotations @@ -30,7 +32,8 @@ class TestPaddingInvariance: # padded-vs-solo AIMNet2 invariant tests/test_model_adapter.py:249-250 # already asserts with atol=1e-2 eV, and ANI2xt's float32 output caps # usable precision at ~float32 ULP (~4e-3 eV) at typical total-energy - # magnitudes per src/Auto3D/batch_opt/ANI2xt_no_rep.py:148-155. 1e-6 would + # magnitudes per the float32-precision Note in ``ANI2xt.forward``'s + # docstring (``src/Auto3D/batch_opt/ANI2xt_no_rep.py``). 1e-6 would # demand sub-ULP reproducibility and flake on a correct model. ANI2x # (torchani, periodic-table indexing) shares ANI2xt's float32 output and # the same -1 species_pad convention, so it gets the same 1e-3 budget diff --git a/tests/test_ranking.py b/tests/test_ranking.py index de7be5f3..795d5caf 100644 --- a/tests/test_ranking.py +++ b/tests/test_ranking.py @@ -47,10 +47,10 @@ def _write_mols_to_sdf(mols: list[Chem.Mol], filepath: str) -> None: class TestConformerRankerWithOptimizedFiltering: - """Tests for ConformerRanker with optimized filtering.""" + """Tests for ConformerRanker's conformer filtering.""" - def test_ranker_with_optimized_filtering_default(self, tmp_path): - """ConformerRanker should use optimized filtering by default.""" + def test_ranker_deduplicates_identical_conformers(self, tmp_path): + """ConformerRanker dedups through the single conformer filter.""" from Auto3D.ranking import ConformerRanker # Create test molecules - all with same SMILES root name. @@ -72,82 +72,66 @@ def test_ranker_with_optimized_filtering_default(self, tmp_path): k=5, ) - # Default should use optimized filtering - assert ranker.use_optimized_filtering is True - results = ranker.run() # With close energies, all in same cluster, should deduplicate to 1 assert len(results) == 1 - def test_ranker_with_legacy_filtering_fallback(self, tmp_path): - """ConformerRanker should support legacy filtering when explicitly requested.""" - from Auto3D.ranking import ConformerRanker + def test_use_optimized_filtering_is_gone_not_silently_ignored(self, tmp_path): + """The flag that selected between two filter implementations is removed. - # Same structure and near-identical energy (within the duplicate energy - # tolerance) -> one unique structure. - mol1 = _create_mol_with_energy("C", -10.0, "mol_1") - mol2 = _create_mol_with_energy("C", -10.005, "mol_2") + There were two conformer filters with the same duplicate criterion, + chosen by ``use_optimized_filtering``, and they had drifted on malformed + input. One filter now, and the flag must raise rather than be swallowed + by ``**kwargs`` -- a caller passing ``False`` deserves to learn the + legacy path is gone instead of silently getting the other one. + """ + from Auto3D.ranking import ConformerRanker input_path = str(tmp_path / "input.sdf") - output_path = str(tmp_path / "output.sdf") - _write_mols_to_sdf([mol1, mol2], input_path) - - ranker = ConformerRanker( - input_path=input_path, - out_path=output_path, - threshold=0.3, - k=5, - use_optimized_filtering=False, - ) - - assert ranker.use_optimized_filtering is False - - results = ranker.run() - # Same behavior, one unique structure - assert len(results) == 1 - - def test_ranker_optimized_vs_legacy_produce_same_results(self, tmp_path): - """Optimized and legacy filtering should produce equivalent results. - - With a large energy_cluster_window, optimized should behave like legacy. + _write_mols_to_sdf([_create_mol_with_energy("C", -10.0, "mol_1")], input_path) + + with pytest.raises(TypeError, match="use_optimized_filtering"): + ConformerRanker( + input_path=input_path, + out_path=str(tmp_path / "output.sdf"), + threshold=0.3, + k=5, + use_optimized_filtering=False, + ) + + def test_a_single_cluster_gives_the_same_answer(self, tmp_path): + """A window wide enough to make one cluster must not change the result. + + That equivalence is the entire justification for partitioning the energy + axis at all, and it used to be asserted by comparing against the legacy + all-pairs filter (deleted in 4.1.0). Asserted directly now. """ from Auto3D.ranking import ConformerRanker # Identical molecules (same structure AND near-identical energy, within # the duplicate energy tolerance) - these should be deduplicated. - mol1 = _create_mol_with_energy("CCCC", -10.0, "a_1") - mol2 = _create_mol_with_energy("CCCC", -10.005, "a_2") # same structure & energy - mol3 = _create_mol_with_energy("CCCC", -10.008, "a_3") # same structure & energy + mols = [ + _create_mol_with_energy("CCCC", -10.0, "a_1"), + _create_mol_with_energy("CCCC", -10.005, "a_2"), + _create_mol_with_energy("CCCC", -10.008, "a_3"), + ] input_path = str(tmp_path / "input.sdf") - output_optimized = str(tmp_path / "output_optimized.sdf") - output_legacy = str(tmp_path / "output_legacy.sdf") - _write_mols_to_sdf([mol1, mol2, mol3], input_path) - - # Test with optimized filtering - use large energy window to match legacy behavior - ranker_optimized = ConformerRanker( - input_path=input_path, - out_path=output_optimized, - threshold=0.3, - k=5, - use_optimized_filtering=True, - energy_cluster_window=100.0, # Large window = single cluster = legacy behavior - ) - results_optimized = ranker_optimized.run() + _write_mols_to_sdf(mols, input_path) - # Test with legacy filtering - ranker_legacy = ConformerRanker( - input_path=input_path, - out_path=output_legacy, - threshold=0.3, - k=5, - use_optimized_filtering=False, - ) - results_legacy = ranker_legacy.run() + counts = [] + for name, window in (("default", None), ("single_cluster", 100.0)): + kwargs = {} if window is None else {"energy_cluster_window": window} + counts.append(len(ConformerRanker( + input_path=input_path, + out_path=str(tmp_path / f"output_{name}.sdf"), + threshold=0.3, + k=5, + **kwargs, + ).run())) - # Should have same number of results - all identical molecules deduplicated to 1 - assert len(results_optimized) == len(results_legacy) - assert len(results_optimized) == 1 # All identical molecules should be deduplicated + assert counts[0] == counts[1] + assert counts[0] == 1 # all identical molecules deduplicate to one def test_energy_cluster_window_parameter(self, tmp_path): """Ranker should accept energy_cluster_window parameter for optimized filtering.""" @@ -198,7 +182,6 @@ def test_top_k_with_optimized_filtering(self, tmp_path): out_path=output_path, threshold=0.3, k=2, - use_optimized_filtering=True, ) df = pd.DataFrame({ @@ -367,7 +350,6 @@ def test_top_window_with_optimized_filtering(self, tmp_path): out_path=output_path, threshold=0.3, window=1.0, # 1 kcal/mol window - use_optimized_filtering=True, ) df = pd.DataFrame({ @@ -789,3 +771,453 @@ def test_enumerate_isomer_false_returns_both_molecules(self, tmp_path): ] assert {m.GetProp("_Name") for m in written} == {key, key_2} assert os.path.getsize(output) > 0 + + +class TestNothingSelectedSaysWhy: + """"No structure converged" used to be the message for every empty group. + + ``ranking`` logged it whether the conformers were dropped for convergence, + for stereochemistry (an optimization that inverted a center, so the geometry + no longer matches the title) or for connectivity (a structure that fell + apart). Two of those three point at the input or the chemistry, and both + were reported as an optimizer convergence problem -- sending the reader to + ``--opt-steps`` and ``--convergence-threshold`` for something neither can + fix. + """ + + @staticmethod + def _group(mols: list[Chem.Mol], name: str = "probe"): + """A one-species ranking group, shaped the way ``run`` builds them.""" + import pandas as pd + + return pd.DataFrame({ + "names": [name] * len(mols), + "energies": [e_tot_ev(m) for m in mols], + "mols": mols, + }) + + @staticmethod + def _ranker(tmp_path, **kwargs): + from Auto3D.ranking import ConformerRanker + + return ConformerRanker( + input_path=str(tmp_path / "in.sdf"), + out_path=str(tmp_path / "out.sdf"), + threshold=0.3, + **kwargs, + ) + + @staticmethod + def _stereo_changed(energy: float, name: str) -> Chem.Mol: + from Auto3D.utils.stereo_check import STEREO_CHANGED_PROP + + mol = _create_mol_with_energy("C/C=C/CCO", energy, name) + mol.SetProp(STEREO_CHANGED_PROP, "true") + return mol + + @staticmethod + def _messages(caplog) -> list[str]: + return [r.getMessage() for r in caplog.records] + + def test_a_stereo_dropped_species_is_not_called_unconverged(self, tmp_path, caplog): + import logging + + mols = [self._stereo_changed(-10.0, "probe_0_0"), + self._stereo_changed(-9.0, "probe_0_1")] + ranker = self._ranker(tmp_path, k=5) + + with caplog.at_level(logging.INFO, logger="Auto3D.ranking"): + assert ranker.top_k(self._group(mols), k=5) == [] + + messages = self._messages(caplog) + assert not any("No structure converged" in m for m in messages), ( + f"conformers dropped for stereochemistry were reported as an " + f"optimizer convergence failure: {messages}" + ) + assert any( + "probe" in m and "stereochemistry" in m for m in messages + ), f"the real reason was never named: {messages}" + + def test_the_k1_fast_path_reports_the_same_reason(self, tmp_path, caplog): + """k=1 bypasses the RMSD dedup entirely; the diagnostic must not + depend on which k the user asked for.""" + import logging + + mols = [self._stereo_changed(-10.0, "probe_0_0")] + ranker = self._ranker(tmp_path, k=1) + + with caplog.at_level(logging.INFO, logger="Auto3D.ranking"): + assert ranker.top_k(self._group(mols), k=1) == [] + + messages = self._messages(caplog) + assert not any("No structure converged" in m for m in messages), messages + assert any("stereochemistry" in m for m in messages), messages + + def test_a_connectivity_dropped_species_names_connectivity(self, tmp_path, caplog): + import logging + + from Auto3D.utils.connectivity import check_connectivity + + def broken(energy: float, name: str) -> Chem.Mol: + mol = _create_mol_with_energy("CC", energy, name) + conf = mol.GetConformer() + pos = conf.GetAtomPosition(0) + conf.SetAtomPosition(0, (pos.x + 5.0, pos.y, pos.z)) + assert check_connectivity(mol) is False, "test premise" + return mol + + mols = [broken(-10.0, "probe_0_0"), broken(-9.0, "probe_0_1")] + ranker = self._ranker(tmp_path, k=5) + + with caplog.at_level(logging.INFO, logger="Auto3D.ranking"): + assert ranker.top_k(self._group(mols), k=5) == [] + + messages = self._messages(caplog) + assert not any("No structure converged" in m for m in messages), messages + assert any("broken or newly formed bonds" in m for m in messages), messages + + def test_the_literal_survives_when_convergence_is_the_sole_reason( + self, tmp_path, caplog + ): + """The inverse, and the reason the assertions above are safe. + + A change that simply stopped emitting "No structure converged" would + satisfy every test above while destroying the message users and their + log-scraping scripts have matched on since Auto3D 1.x. When convergence + IS the sole reason, the exact wording must still appear. + """ + import logging + + mols = [ + _create_mol_with_energy("CCO", -10.0, "probe_0_0", converged=False), + _create_mol_with_energy("CCO", -9.0, "probe_0_1", converged=False), + ] + ranker = self._ranker(tmp_path, k=5) + + with caplog.at_level(logging.INFO, logger="Auto3D.ranking"): + assert ranker.top_k(self._group(mols), k=5) == [] + + assert any( + m == "No structure converged for probe." for m in self._messages(caplog) + ), self._messages(caplog) + + def test_the_literal_is_not_emitted_alongside_another_reason( + self, tmp_path, caplog + ): + """"Only when that is the sole reason": a mixed group must not claim + convergence.""" + import logging + + mols = [ + _create_mol_with_energy("C/C=C/CCO", -10.0, "probe_0_0", converged=False), + self._stereo_changed(-9.0, "probe_0_1"), + ] + ranker = self._ranker(tmp_path, k=5) + + with caplog.at_level(logging.INFO, logger="Auto3D.ranking"): + assert ranker.top_k(self._group(mols), k=5) == [] + + messages = self._messages(caplog) + assert not any("No structure converged" in m for m in messages), messages + # Both reasons are named, so the reader sees the whole accounting. + assert any( + "Converged=false" in m and "stereochemistry" in m for m in messages + ), messages + + def test_a_successful_selection_logs_no_complaint(self, tmp_path, caplog): + """The other inverse: a group that DOES select must stay silent. + + A message emitted unconditionally would pass the "names the reason" + tests above and spam every ordinary run. + """ + import logging + + mols = [_create_mol_with_energy("CCO", -10.0, "probe_0_0")] + ranker = self._ranker(tmp_path, k=5) + + with caplog.at_level(logging.INFO, logger="Auto3D.ranking"): + assert len(ranker.top_k(self._group(mols), k=5)) == 1 + + messages = self._messages(caplog) + assert not any("No structure" in m for m in messages), messages + + def test_top_window_merges_the_window_into_the_same_accounting(self, tmp_path): + """The energy window is the one drop reason ``top_window`` owns. + + It goes into the same run-level tally as the filter's own counts, so + ``run``'s summary is one accounting of the whole selection rather than + two partial ones. + """ + ranker = self._ranker(tmp_path, window=1.0) + + # Two distinct compounds, 5 eV apart: far outside a 1 kcal/mol window. + mols = [ + _create_mol_with_energy("CCO", -10.0, "probe_0_0"), + _create_mol_with_energy("CCCCCCO", -5.0, "probe_0_1"), + ] + kept = ranker.top_window(self._group(mols), window=1.0) + + assert len(kept) == 1, "the second conformer is outside the window" + assert ranker._drop_totals == {"energy_window": 1} + + def test_a_wide_window_records_no_window_drop(self, tmp_path): + """The inverse: a window nothing falls outside of must not report one.""" + ranker = self._ranker(tmp_path, window=1000.0) + mols = [ + _create_mol_with_energy("CCO", -10.0, "probe_0_0"), + _create_mol_with_energy("CCCCCCO", -5.0, "probe_0_1"), + ] + assert len(ranker.top_window(self._group(mols), window=1000.0)) == 2 + assert ranker._drop_totals == {} + + def test_the_run_level_warning_names_the_reasons_that_fired( + self, tmp_path, caplog + ): + """``run``'s "selected 0 structures" warning used to list every reason + it MIGHT have been. + + The text was "N record(s) are marked Converged=false and the rest were + dropped by the connectivity, stereochemistry or energy-window filters" + -- a hand-maintained disjunction that left the reader to work out which + of the three actually happened, on the one message that reaches a direct + API caller's stderr. Here nothing is unconverged and everything is + stereo-changed, so naming convergence would be wrong. + """ + import logging + + from Auto3D.ranking import ConformerRanker + + mols = [self._stereo_changed(-10.0, "probe_0_0"), + self._stereo_changed(-9.0, "probe_0_1")] + input_path = str(tmp_path / "in.sdf") + _write_mols_to_sdf(mols, input_path) + + ranker = ConformerRanker( + input_path=input_path, + out_path=str(tmp_path / "out.sdf"), + threshold=0.3, + k=5, + ) + with caplog.at_level(logging.WARNING, logger="Auto3D.ranking"): + assert ranker.run() == [] + + warnings = [ + r.getMessage() for r in caplog.records if r.levelno >= logging.WARNING + ] + assert any("Selected 0 structures from 2 record(s)" in m for m in warnings), ( + warnings + ) + assert any( + "changed stereochemistry during optimization" in m for m in warnings + ), warnings + assert not any("Converged=false" in m for m in warnings), ( + f"nothing was unconverged, yet convergence was named: {warnings}" + ) + + def test_the_run_level_warning_still_names_convergence_when_it_applies( + self, tmp_path, caplog + ): + """The inverse: records dropped before grouping (for convergence, or + because RDKit could not parse them) are counted in the same tally.""" + import logging + + from Auto3D.ranking import ConformerRanker + + mols = [ + _create_mol_with_energy("CCO", -10.0, "probe_0_0", converged=False), + _create_mol_with_energy("CCO", -9.0, "probe_0_1", converged=False), + ] + input_path = str(tmp_path / "in.sdf") + _write_mols_to_sdf(mols, input_path) + + ranker = ConformerRanker( + input_path=input_path, + out_path=str(tmp_path / "out.sdf"), + threshold=0.3, + k=5, + ) + with caplog.at_level(logging.WARNING, logger="Auto3D.ranking"): + assert ranker.run() == [] + + warnings = [ + r.getMessage() for r in caplog.records if r.levelno >= logging.WARNING + ] + assert any("2 marked Converged=false" in m for m in warnings), warnings + + def test_the_tally_is_reset_between_runs(self, tmp_path): + """A ranker reused for a second file must not report the first's drops.""" + from Auto3D.ranking import ConformerRanker + + # Stereo-changed records pass run()'s own convergence check and are + # dropped by the FILTER, so they land in the tally -- unlike an + # unconverged record, which run() drops before grouping and counts + # separately. Getting that wrong makes this test vacuous. + dirty = str(tmp_path / "dirty.sdf") + _write_mols_to_sdf( + [self._stereo_changed(-10.0, "a_0_0"), + self._stereo_changed(-9.0, "a_0_1")], + dirty, + ) + clean = str(tmp_path / "clean.sdf") + _write_mols_to_sdf([_create_mol_with_energy("CCO", -10.0, "b_0_0")], clean) + + ranker = ConformerRanker( + input_path=dirty, out_path=str(tmp_path / "out.sdf"), threshold=0.3, k=5, + ) + assert ranker.run() == [] + assert ranker._drop_totals == {"stereochemistry": 2}, "test premise" + + ranker.input_path = clean + assert len(ranker.run()) == 1 + assert ranker._drop_totals == {}, ( + "the second run reported the first run's drops" + ) + + +class TestSelectorDispatchRegistry: + """``run`` dispatches through a registry checked against ``config.py``. + + It used to be a hand-written ``if self.k: ... elif self.window: ...``. A + third selector added to ``Auto3D.config.SELECTOR_FIELDS`` would then be + accepted by ``Auto3DOptions``, accepted by ``CLIConfig``, accepted by + ``check_selectors_mutually_exclusive`` -- and silently ignored here, falling + through to "Parameter k or window needs to be specified" even though the + user had specified one. + """ + + def test_the_registry_matches_the_authoritative_field_list(self): + from Auto3D.config import SELECTOR_FIELDS + from Auto3D.ranking import _SELECTORS + + assert set(_SELECTORS) == set(SELECTOR_FIELDS) + + def test_every_mapped_method_exists(self): + from Auto3D.ranking import _SELECTORS, ConformerRanker + + for field, method in _SELECTORS.items(): + assert callable(getattr(ConformerRanker, method, None)), field + + def test_a_selector_in_config_with_nothing_wired_to_it_is_refused(self): + """The point of the check: adding a field to config without wiring a + method here must be impossible to miss. + + This is the exact call ``Auto3D.ranking`` makes at import, with the + field list a developer would have just extended. + """ + from Auto3D.ranking import ( + _SELECTORS, + ConformerRanker, + _verify_selector_registry, + ) + + with pytest.raises(ImportError, match="out of step with"): + _verify_selector_registry( + _SELECTORS, ("k", "window", "percentile"), ConformerRanker + ) + + def test_a_registry_entry_this_module_does_not_implement_is_refused(self): + """A typo in a method name passes the set comparison, then raises + AttributeError from inside ``run`` -- after the whole input has been + read and grouped.""" + from Auto3D.ranking import ConformerRanker, _verify_selector_registry + + with pytest.raises(ImportError, match="not a method of ConformerRanker"): + _verify_selector_registry( + {"k": "top_k", "window": "top_windwo"}, + ("k", "window"), + ConformerRanker, + ) + + def test_the_real_registry_passes_its_own_check(self): + """The inverse: a check that refused everything would satisfy both + tests above and make ``import Auto3D.ranking`` impossible -- so assert + the shipped configuration is accepted.""" + from Auto3D.config import SELECTOR_FIELDS + from Auto3D.ranking import ( + _SELECTORS, + ConformerRanker, + _verify_selector_registry, + ) + + _verify_selector_registry(_SELECTORS, SELECTOR_FIELDS, ConformerRanker) + + def test_a_third_selector_is_dispatched_not_ignored(self, tmp_path, monkeypatch): + """The defect the registry exists to prevent, exercised end to end. + + With the old hand-written ``if self.k: ... elif self.window: ...``, a + selector added to ``SELECTOR_FIELDS`` and wired here was still ignored: + ``run`` consulted only the two names baked into the chain, so a user who + specified the new selector got "Parameter k or window needs to be + specified" for a parameter they had specified. Reverting ``run`` to that + chain must fail this test. + """ + import Auto3D.ranking as ranking + from Auto3D.ranking import ConformerRanker + + calls = [] + + def top_percentile(self, df_group, percentile): + calls.append(percentile) + selected = list(df_group["mols"])[:1] + for mol in selected: + # Every real selector sets this; run() converts it on the way out. + mol.SetProp("E_rel(eV)", "0.0") + return selected + + monkeypatch.setattr( + ranking, "SELECTOR_FIELDS", ("k", "window", "percentile"), raising=True + ) + monkeypatch.setattr( + ranking, + "_SELECTORS", + {"k": "top_k", "window": "top_window", "percentile": "top_percentile"}, + ) + monkeypatch.setattr( + ConformerRanker, "top_percentile", top_percentile, raising=False + ) + + mols = [_create_mol_with_energy("CCO", -10.0, "probe_0_0")] + input_path = str(tmp_path / "in.sdf") + _write_mols_to_sdf(mols, input_path) + + ranker = ConformerRanker( + input_path=input_path, + out_path=str(tmp_path / "out.sdf"), + threshold=0.3, + ) + ranker.percentile = 90.0 + + results = ranker.run() + + assert calls == [90.0], ( + "the third selector was never dispatched to; run() consulted a " + "hard-coded list of selector names instead of the registry" + ) + assert len(results) == 1 + + def test_k_routes_to_top_k_and_window_to_top_window(self, tmp_path): + """The registry must actually be what dispatch consults.""" + from Auto3D.ranking import ConformerRanker + + mols = [_create_mol_with_energy("CCO", -10.0, "probe_0_0")] + input_path = str(tmp_path / "in.sdf") + _write_mols_to_sdf(mols, input_path) + + for field, expected, value in (("k", "top_k", 1), ("window", "top_window", 5.0)): + called = [] + ranker = ConformerRanker( + input_path=input_path, + out_path=str(tmp_path / f"out_{field}.sdf"), + threshold=0.3, + **{field: value}, + ) + for method in ("top_k", "top_window"): + original = getattr(ranker, method) + + def spy(*args, _m=method, _o=original, _log=called, **kwargs): + _log.append(_m) + return _o(*args, **kwargs) + + setattr(ranker, method, spy) + ranker.run() + assert called == [expected], f"{field} dispatched to {called}" diff --git a/tests/test_stereo_postopt.py b/tests/test_stereo_postopt.py index 1bee6be3..8ed2d8df 100644 --- a/tests/test_stereo_postopt.py +++ b/tests/test_stereo_postopt.py @@ -14,7 +14,6 @@ from Auto3D.filtering import filter_unique_optimized from Auto3D.ranking import ConformerRanker -from Auto3D.filtering import filter_unique from Auto3D.utils.stereo_check import ( STEREO_CHANGED_PROP, apply_optimized_coords, @@ -176,12 +175,22 @@ def test_filter_unique_optimized_drops_the_changed_record(self): assert len(result) == 1, f"expected only the preserved record: {len(result)}" assert result[0].GetProp("E_tot") == "-1.0" - def test_filter_unique_drops_the_changed_record(self): - kept = _optimized(-1.0, changed=False) - dropped = _optimized(-2.0, changed=True) - result = filter_unique([dropped, kept], crit=0.3) - assert len(result) == 1, f"expected only the preserved record: {len(result)}" - assert result[0].GetProp("E_tot") == "-1.0" + def test_the_filter_reports_stereochemistry_as_the_drop_reason(self): + """Not just that it dropped, but that it says why. + + Until 4.1.0 the two conformer filters returned a bare list, so + ``ranking`` reported a stereo-changed species as "No structure + converged" -- pointing the reader at the optimizer settings for a + problem in the input's stereo definitions. + """ + from Auto3D.filtering import filter_conformers + + result = filter_conformers( + [_optimized(-2.0, changed=True), _optimized(-1.0, changed=False)], + rmsd_threshold=0.3, + ) + assert result.dropped == {"stereochemistry": 1} + assert result.reasons == ("stereochemistry",) def test_top_k_one_skips_the_changed_lowest_energy_record(self): """k=1 takes a fast path that bypasses the RMSD filters entirely.""" @@ -207,4 +216,3 @@ def test_unmarked_records_still_survive_every_filter(self): """No regression for molecules that never went through the check.""" mols = [_optimized(-1.0, changed=None), _optimized(-2.0, changed=None)] assert len(filter_unique_optimized(mols, rmsd_threshold=0.3)) == 2 - assert len(filter_unique(mols, crit=0.3)) == 2 diff --git a/tests/test_utils_stereochemistry.py b/tests/test_utils_stereochemistry.py index 0a66bfc5..01c9dea9 100644 --- a/tests/test_utils_stereochemistry.py +++ b/tests/test_utils_stereochemistry.py @@ -9,6 +9,7 @@ import pytest from rdkit import Chem +from Auto3D.exceptions import InputValidationError from Auto3D.utils.stereochemistry import ( amend_configuration, amend_configuration_w, @@ -202,6 +203,30 @@ def test_multiple_centers_all_inverted(self): for orig, res in zip(orig_info, result_info): assert orig != res + def test_three_centers_all_inverted(self): + """M60 regression: 3+ centers used to be handled by a loop reading a + variable (key2) set inside the loop's else-branch and read again + after the loop -- correct only by Python's lack of block scoping. + A single-pass rewrite must still invert every center in order.""" + smi = "C[C@H](O)[C@@H](F)[C@H](Cl)Br" + result = create_enantiomer(smi) + assert result == "C[C@@H](O)[C@H](F)[C@@H](Cl)Br" + orig_info = list(get_stereo_info(smi).values()) + result_info = list(get_stereo_info(result).values()) + assert len(result_info) == len(orig_info) == 3 + for orig, res in zip(orig_info, result_info): + assert orig != res + + def test_four_centers_all_inverted(self): + """Same as above, one more center, to rule out an off-by-one at the + boundary between the removed len(keys)==1 special case and the + general multi-key path.""" + smi = "C[C@H](O)[C@@H](F)[C@H](Cl)[C@@H](Br)I" + result = create_enantiomer(smi) + assert result == "C[C@@H](O)[C@H](F)[C@@H](Cl)[C@H](Br)I" + result_info = list(get_stereo_info(result).values()) + assert len(result_info) == 4 + class TestCheckValue: """Tests for the check_value() function.""" @@ -281,6 +306,32 @@ def test_remove_enantiomers_no_stereo(self): os.unlink(inpath) os.unlink(outpath) + def test_remove_enantiomers_tolerates_blank_lines(self, tmp_path): + """M59: switched to iter_smi_records, which skips blank/comment lines. + + The old hand-rolled parser did `vals = line.split()` then indexed + vals[0]/vals[1] unconditionally, so a blank line raised a bare + IndexError and aborted the whole function. + """ + inpath = tmp_path / "in.smi" + inpath.write_text("CCO mol1_1\n\n# a comment\nCCCO mol2_1\n") + outpath = tmp_path / "out.smi" + + result = remove_enantiomers(str(inpath), str(outpath)) + + assert "mol1" in result + assert "mol2" in result + + def test_remove_enantiomers_rejects_missing_id(self, tmp_path): + """A non-blank, non-comment line with no ID must still fail loudly, + just as an InputValidationError rather than a bare IndexError.""" + inpath = tmp_path / "in.smi" + inpath.write_text("CCO\n") + outpath = tmp_path / "out.smi" + + with pytest.raises(InputValidationError): + remove_enantiomers(str(inpath), str(outpath)) + class TestAmendConfiguration: """Tests for the amend_configuration() function.""" @@ -317,6 +368,39 @@ def test_amend_configuration_w_writes_file(self): finally: os.unlink(path) + def test_amend_configuration_tolerates_blank_and_comment_lines(self, tmp_path): + """M59: switched to iter_smi_records. The old parser did + `tuple(line.strip().split())`, which raises "not enough values to + unpack" on a blank line -- aborting the whole function -- and "too + many values to unpack" on a line with a third whitespace column. + """ + path = tmp_path / "in.smi" + path.write_text( + "C[C@H](O)F mol_1\n\n# comment\nC[C@@H](O)F mol_2\n" + ) + + result = amend_configuration(str(path)) + + assert "mol" in result + + def test_amend_configuration_tolerates_extra_column(self, tmp_path): + """A trailing whitespace-separated column beyond SMILES+ID must not + raise, matching every other consumer of this format.""" + path = tmp_path / "in.smi" + path.write_text("C[C@H](O)F mol_1 extra_column\n") + + result = amend_configuration(str(path)) + + assert "mol" in result + + def test_amend_configuration_rejects_missing_id(self, tmp_path): + """A non-blank, non-comment line with no ID must still fail loudly.""" + path = tmp_path / "in.smi" + path.write_text("C[C@H](O)F\n") + + with pytest.raises(InputValidationError): + amend_configuration(str(path)) + class TestIntegration: """Integration tests using real molecular examples.""" diff --git a/tests/test_utils_validation.py b/tests/test_utils_validation.py index a1643907..67d86e15 100644 --- a/tests/test_utils_validation.py +++ b/tests/test_utils_validation.py @@ -6,7 +6,6 @@ from rdkit import Chem from Auto3D.config import Auto3DOptions -from Auto3D.filtering import filter_unique from Auto3D.utils.connectivity import check_connectivity from Auto3D.utils.validation import ( check_input, @@ -150,113 +149,6 @@ def test_check_connectivity_example_sdf(self): assert check_connectivity(mol) is True -class TestFilterUnique: - """Tests for filter_unique function.""" - - def test_filter_unique_removes_unconverged(self): - """Test that filter_unique removes unconverged structures.""" - supp = Chem.SDMolSupplier(path_example_sdf, removeHs=False) - mols = [mol for mol in supp if mol is not None] - - # Set all mols as unconverged - for mol in mols: - mol.SetProp("Converged", "False") - - result = filter_unique(mols) - assert len(result) == 0 - - def test_filter_unique_keeps_converged(self): - """Test that filter_unique keeps converged structures.""" - supp = Chem.SDMolSupplier(path_example_sdf, removeHs=False) - mols = [mol for mol in supp if mol is not None] - - # Set all mols as converged - for mol in mols: - mol.SetProp("Converged", "True") - - result = filter_unique(mols) - # Should keep at least one unique structure - assert len(result) >= 1 - - def test_filter_unique_removes_duplicates(self): - """Test that filter_unique removes similar structures.""" - supp = Chem.SDMolSupplier(path_example_sdf, removeHs=False) - mols = [mol for mol in supp if mol is not None] - - if len(mols) > 0: - # Duplicate the first molecule - mol = mols[0] - mol.SetProp("Converged", "True") - duplicate = Chem.RWMol(mol) - duplicate.SetProp("Converged", "True") - - result = filter_unique([mol, duplicate], crit=0.3) - # Only one should remain - assert len(result) == 1 - - def test_filter_unique_keeps_records_with_no_converged_property(self): - """Absence of the property is not a failed optimization. - - Only ``batchopt`` writes ``Converged``; an ``opt_geometry`` output, an - ORCA/Gaussian export or a hand-built conformer set carries none. - Treating that as "did not converge" deleted every record. The - consequence asserted here is that such a file filters exactly the same - as one whose records all say Converged=True. - """ - supp = Chem.SDMolSupplier(path_example_sdf, removeHs=False) - flagged = [mol for mol in supp if mol is not None] - for mol in flagged: - mol.SetProp("Converged", "True") - expected = len(filter_unique(flagged)) - assert expected >= 1, "test premise: the flagged file must keep something" - - supp = Chem.SDMolSupplier(path_example_sdf, removeHs=False) - unflagged = [mol for mol in supp if mol is not None] - for mol in unflagged: - mol.ClearProp("Converged") - assert not mol.HasProp("Converged") - - result = filter_unique(unflagged) - assert len(result) == expected, ( - f"{len(unflagged)} record(s) with no 'Converged' property kept " - f"{len(result)}, but the same records marked Converged=True keep " - f"{expected}" - ) - - def test_filter_unique_custom_threshold(self): - """A tighter RMSD threshold must keep MORE structures than a looser - one -- not merely "at least as many", which passes even when the two - thresholds produce identical results (as the fixture in - ``path_example_sdf`` does: two molecules of different sizes, so - ``species_key`` alone already keeps both regardless of ``crit``). - - Constructs two conformers of ONE molecule whose RMSD sits strictly - between the two thresholds, so equality cannot pass silently. - """ - mol1 = Chem.AddHs(Chem.MolFromSmiles("CCCCCCCC")) # octane: flexible - mol2 = Chem.Mol(mol1) - from rdkit.Chem import AllChem, rdMolAlign - AllChem.EmbedMolecule(mol1, randomSeed=1) - AllChem.EmbedMolecule(mol2, randomSeed=99) - AllChem.MMFFOptimizeMolecule(mol1) - AllChem.MMFFOptimizeMolecule(mol2) - mol1.SetProp("Converged", "True") - mol2.SetProp("Converged", "True") - - rmsd = rdMolAlign.GetBestRMS(Chem.RemoveHs(mol1), Chem.RemoveHs(mol2)) - assert rmsd > 0.05, "test premise: conformers must be geometrically distinct" - - crit_strict = rmsd / 2 # below the actual RMSD -> kept separate - crit_lenient = rmsd * 2 # above the actual RMSD -> merged - - result_strict = filter_unique([mol1, mol2], crit=crit_strict) - result_lenient = filter_unique([mol1, mol2], crit=crit_lenient) - - assert len(result_strict) == 2, "strict threshold must not merge distinct conformers" - assert len(result_lenient) == 1, "lenient threshold must merge near-identical conformers" - assert len(result_strict) > len(result_lenient) - - def _options(**overrides): """Build the ``Auto3DOptions`` ``check_valid_configuration`` now takes.