diff --git a/docs/unified-api.md b/docs/unified-api.md index 3b6aff7..4500fea 100644 --- a/docs/unified-api.md +++ b/docs/unified-api.md @@ -94,6 +94,15 @@ Paths in CLI specifications are relative to the specification file. Model names remain identifiers unless they resolve to local checkpoint files. Python paths are relative to the caller's working directory. +For Zeo++ surface area, volume, and PSD, `channel_radius` controls accessibility +and `probe_radius` controls sampling. The probe radius must not exceed the +channel radius; equal radii remain the default. Both APIs pass these values in +Zeo++'s documented channel/probe order. All requested outputs must be present, +complete, and finite. Collection accepts multiline channel output and the +default `.psd_histo` histogram filename, retaining explicitly named legacy +`.psd` files and existing result keys. See the +[Zeo++ command and output examples](https://www.zeoplusplus.org/examples.html). + ```bash matkit pores --spec pores.json --outdir runs/pores matkit prepare --spec pores.json --outdir runs/prepared diff --git a/src/matkit/api/adapters.py b/src/matkit/api/adapters.py index 873d9fb..07c023a 100644 --- a/src/matkit/api/adapters.py +++ b/src/matkit/api/adapters.py @@ -289,47 +289,30 @@ def external_command(root, execution, engine, arguments): def run_pores(root, request, execution): - args = ["-ha"] if request.high_accuracy else [] - args.extend(["-r", request.radii_file]) - for analysis in request.analyses: - args.append(f"-{analysis}") - if analysis in {"sa", "vol", "psd"}: - args.extend( - [ - str(request.probe_radius), - str(request.channel_radius), - str(request.num_samples), - ] - ) - elif analysis == "chan": - args.append(str(request.probe_radius)) - args.append("structure.cif") + from matkit.zeopp.zeopp import _analysis_arguments + + args = _analysis_arguments( + request.analyses, + request.probe_radius, + request.channel_radius, + request.num_samples, + request.high_accuracy, + request.radii_file, + "structure.cif", + ) external_command(root, execution, "zeopp", args) return parse_pores(root, request) def parse_pores(root, request): - from matkit.zeopp.zeopp import _PARSERS - - required = { - "res": {"Di", "Df", "Dif"}, - "sa": {"ASA", "NASA", "density", "unitcell_volume"}, - "vol": {"AV", "NAV", "density", "unitcell_volume"}, - "psd": {"bin_lower", "counts"}, - "chan": {"num_channels", "dimensionalities"}, + from matkit.zeopp.zeopp import _output_path, _parse_output + + results = { + analysis: _parse_output( + _output_path(root / "work", "structure", analysis), analysis + ) + for analysis in request.analyses } - results = {} - for analysis in request.analyses: - data = _PARSERS[analysis](root / "work" / f"structure.{analysis}") - if not required[analysis] <= data.keys(): - raise ValueError(f"Incomplete requested Zeo++ analysis: {analysis}") - if analysis == "psd" and not data["counts"]: - raise ValueError("Empty pore-size distribution") - if analysis == "chan" and data["num_channels"] != len( - data["dimensionalities"] - ): - raise ValueError("Channel count does not match dimensionalities") - results[analysis] = data return ( PorePayload(results=results), [ diff --git a/src/matkit/api/models.py b/src/matkit/api/models.py index d147a38..a8fc867 100644 --- a/src/matkit/api/models.py +++ b/src/matkit/api/models.py @@ -231,6 +231,10 @@ class PoreRequest(RequestBase): def unique_analyses(self): if len(set(self.analyses)) != len(self.analyses): raise ValueError("Analyses must be unique") + if {"sa", "vol", "psd"}.intersection( + self.analyses + ) and self.probe_radius > self.channel_radius: + raise ValueError("probe_radius must not exceed channel_radius") return self diff --git a/src/matkit/zeopp/zeopp.py b/src/matkit/zeopp/zeopp.py index 592da41..7b3681b 100644 --- a/src/matkit/zeopp/zeopp.py +++ b/src/matkit/zeopp/zeopp.py @@ -3,6 +3,8 @@ import logging from concurrent.futures import ThreadPoolExecutor import json +import math +from numbers import Integral, Real from pathlib import Path import shutil import subprocess @@ -14,6 +16,61 @@ VALID_ANALYSES = {"res", "sa", "vol", "psd", "chan"} +_OUTPUT_SUFFIXES = {analysis: (f".{analysis}",) for analysis in VALID_ANALYSES} +# Zeo++ uses this suffix by default; retain explicitly named legacy files. +_OUTPUT_SUFFIXES["psd"] = (".psd_histo", ".psd") + + +def _analysis_arguments( + analyses, + probe_radius, + channel_radius, + num_samples, + high_accuracy, + radii_file, + structure_file, +): + """Build the engine syntax shared by legacy and unified execution.""" + if not analyses or set(analyses) - VALID_ANALYSES: + raise ValueError(f"Invalid analysis types: {analyses}") + if len(set(analyses)) != len(analyses): + raise ValueError("Analyses must be unique") + for name, value in ( + ("probe_radius", probe_radius), + ("channel_radius", channel_radius), + ): + if ( + isinstance(value, bool) + or not isinstance(value, Real) + or not math.isfinite(value) + or value <= 0 + ): + raise ValueError(f"{name} must be positive and finite") + if ( + isinstance(num_samples, bool) + or not isinstance(num_samples, Integral) + or num_samples < 1 + ): + raise ValueError("num_samples must be a positive integer") + sampled = {"sa", "vol", "psd"} + if sampled.intersection(analyses) and probe_radius > channel_radius: + raise ValueError("probe_radius must not exceed channel_radius") + args = ["-ha"] if high_accuracy else [] + args.extend(["-r", str(radii_file)]) + for analysis in analyses: + args.append(f"-{analysis}") + if analysis in sampled: + args.extend(map(str, (channel_radius, probe_radius, num_samples))) + elif analysis == "chan": + args.append(str(probe_radius)) + return [*args, str(structure_file)] + + +def _output_path(directory: Path, stem: str, analysis: str) -> Path: + candidates = [ + directory / f"{stem}{suffix}" for suffix in _OUTPUT_SUFFIXES[analysis] + ] + return next((path for path in candidates if path.is_file()), candidates[0]) def _find_network_binary(network_path: str | None = None) -> str: @@ -181,16 +238,34 @@ def _parse_psd(filepath: Path) -> dict: counts = [] for line in lines: parts = line.strip().split() - if len(parts) >= 2: - try: - low = float(parts[0]) - count = float(parts[1]) - bin_lower.append(low) - counts.append(count) - except ValueError: - continue + if not parts: + continue + try: + low = float(parts[0]) + except ValueError: + continue # Zeo++ histogram headers are not numeric rows. + if len(parts) < 2: + raise ValueError("Incomplete PSD histogram row") + values = [float(value) for value in parts] + if not all(math.isfinite(value) for value in values): + raise ValueError("PSD histogram values must be finite") + count = values[1] + if low < 0 or count < 0: + raise ValueError("PSD bins and counts must be nonnegative") + if bin_lower and low <= bin_lower[-1]: + raise ValueError("PSD bins must be strictly increasing") + bin_lower.append(low) + counts.append(count) + + if not counts: + raise ValueError("Empty pore-size distribution") bin_size = bin_lower[1] - bin_lower[0] if len(bin_lower) > 1 else 0.0 + if any( + not math.isclose(b - a, bin_size, rel_tol=1e-6, abs_tol=1e-9) + for a, b in zip(bin_lower, bin_lower[1:]) + ): + raise ValueError("PSD histogram requires uniformly spaced bins") bin_upper = [b + bin_size for b in bin_lower] return { @@ -215,12 +290,41 @@ def _parse_chan(filepath: Path) -> dict: Dict with num_channels and dimensionalities list. """ try: - text = filepath.read_text().strip() - parts = text.split() + lines = [ + line.strip() + for line in filepath.read_text().splitlines() + if line.strip() and not line.lstrip().startswith("#") + ] + parts = lines[0].split() channels_idx = parts.index("channels") num_channels = int(parts[channels_idx - 1]) dim_idx = parts.index("dimensionality") dimensionalities = [int(d) for d in parts[dim_idx + 1 :]] + if ( + num_channels < 0 + or num_channels != len(dimensionalities) + or any(d not in {1, 2, 3} for d in dimensionalities) + or parts[channels_idx + 1 : dim_idx] != ["identified", "of"] + ): + raise ValueError("Invalid channel count or dimensionalities") + channels = set() + for line in lines[1:]: + row = line.split() + if row[0] != "Channel": + continue # Optional summary rows do not list dimensionalities. + if len(row) != 5: + raise ValueError("Incomplete channel detail row") + index = int(row[1]) + diameters = [float(value) for value in row[2:]] + if index in channels or not 0 <= index < num_channels: + raise ValueError("Invalid or duplicate channel index") + if any(not math.isfinite(d) or d < 0 for d in diameters): + raise ValueError( + "Channel diameters must be finite and nonnegative" + ) + channels.add(index) + if channels and len(channels) != num_channels: + raise ValueError("Incomplete channel detail rows") return { "num_channels": num_channels, "dimensionalities": dimensionalities, @@ -238,6 +342,25 @@ def _parse_chan(filepath: Path) -> dict: } +def _parse_output(path: Path, analysis: str) -> dict: + """Require complete, finite results in every public Zeo++ interface.""" + required = { + "res": {"Di", "Df", "Dif"}, + "sa": {"ASA", "NASA", "density", "unitcell_volume"}, + "vol": {"AV", "NAV", "density", "unitcell_volume"}, + "psd": {"bin_lower", "counts"}, + "chan": {"num_channels", "dimensionalities"}, + } + data = _PARSERS[analysis](path) + if not required[analysis] <= data.keys(): + raise ValueError(f"Incomplete requested Zeo++ analysis: {analysis}") + try: + json.dumps(data, allow_nan=False) + except ValueError as exc: + raise ValueError(f"Nonfinite Zeo++ analysis: {analysis}") from exc + return data + + def get_output_data( output_path: str, analyses: list[str] | None = None, @@ -247,7 +370,7 @@ def get_output_data( Args: output_path: Path to directory containing Zeo++ output files. analyses: Which analyses to parse. If None, auto-detects from - available files (.res, .sa, .vol, .psd, .chan). + available files (.res, .sa, .vol, .psd_histo, .psd, .chan). Returns: Dict with 'success' key and per-analysis result sub-dicts. @@ -271,17 +394,27 @@ def get_output_data( if outdir.is_file(): ext = outdir.suffix.lstrip(".") + if ext == "psd_histo": + ext = "psd" + if analyses is not None and set(analyses) != {ext}: + raise ValueError("File does not contain all requested analyses") if ext in _PARSERS: - results[ext] = _PARSERS[ext](outdir) + results[ext] = _parse_output(outdir, ext) results["success"] = True return results # Directory mode: find output files detect = analyses if analyses is not None else list(VALID_ANALYSES) for analysis in detect: - matches = list(outdir.glob(f"*.{analysis}")) + matches = [ + path + for suffix in _OUTPUT_SUFFIXES[analysis] + for path in sorted(outdir.glob(f"*{suffix}")) + ] if matches: - results[analysis] = _PARSERS[analysis](matches[0]) + results[analysis] = _parse_output(matches[0], analysis) + elif analyses is not None: + raise ValueError(f"Missing requested Zeo++ analysis: {analysis}") if len(results) > 1: # has at least one analysis result beyond 'success' results["success"] = True @@ -367,47 +500,18 @@ def run_zeopp( rad_dest = workdir / radii_path.name shutil.copy(radii_path, rad_dest) - # Build command - cmd = [binary] - - if ha: - cmd.append("-ha") - - cmd.extend(["-r", str(rad_dest)]) - - for analysis in analyses: - if analysis == "res": - cmd.extend(["-res"]) - elif analysis == "sa": - cmd.extend( - [ - "-sa", - str(probe_radius), - str(chan_radius), - str(num_samples), - ] - ) - elif analysis == "vol": - cmd.extend( - [ - "-vol", - str(probe_radius), - str(chan_radius), - str(num_samples), - ] - ) - elif analysis == "psd": - cmd.extend( - [ - "-psd", - str(probe_radius), - str(chan_radius), - str(num_samples), - ] - ) - elif analysis == "chan": - cmd.extend(["-chan", str(probe_radius)]) - cmd.append(str(cif_dest)) + cmd = [ + binary, + *_analysis_arguments( + analyses, + probe_radius, + chan_radius, + num_samples, + ha, + rad_dest, + cif_dest, + ), + ] proc = subprocess.run(cmd, capture_output=True, text=True) if proc.returncode != 0: @@ -420,9 +524,12 @@ def run_zeopp( result = {"success": False, "results": {}, "error": None} stem = cif_dest.stem for analysis in analyses: - out_file = workdir / f"{stem}.{analysis}" - if out_file.exists(): - result["results"][analysis] = _PARSERS[analysis](out_file) + out_file = _output_path(workdir, stem, analysis) + if not out_file.is_file(): + raise ValueError( + f"Missing requested Zeo++ analysis: {analysis}" + ) + result["results"][analysis] = _parse_output(out_file, analysis) result["success"] = len(result["results"]) > 0 return result diff --git a/tests/data/zeopp/test_structure.chan b/tests/data/zeopp/test_structure.chan index 3a1ad0f..02886ef 100644 --- a/tests/data/zeopp/test_structure.chan +++ b/tests/data/zeopp/test_structure.chan @@ -1 +1,3 @@ test_structure 2 channels identified of dimensionality 3 3 +Channel 0 18.569 8.023 10.928 +Channel 1 12.0 6.0 8.0 diff --git a/tests/fixtures/fake_engine.py b/tests/fixtures/fake_engine.py index b4860b0..5a01135 100644 --- a/tests/fixtures/fake_engine.py +++ b/tests/fixtures/fake_engine.py @@ -17,7 +17,8 @@ if f"-{analysis}" in sys.argv and not ( "--partial" in sys.argv and analysis == "sa" ): - Path(f"structure.{analysis}").write_bytes( + suffix = "psd_histo" if analysis == "psd" else analysis + Path(f"structure.{suffix}").write_bytes( (fixtures / f"test_structure.{analysis}").read_bytes() ) elif mode == "graspa": diff --git a/tests/test_api_zeopp.py b/tests/test_api_zeopp.py new file mode 100644 index 0000000..41d8dc2 --- /dev/null +++ b/tests/test_api_zeopp.py @@ -0,0 +1,176 @@ +"""Engine syntax and documented output shapes, without an installed Zeo++.""" + +from pathlib import Path +import shutil +from types import SimpleNamespace + +import pytest + +from matkit.api import PoreRequest, StructureRef, run +from matkit.api import adapters +from matkit.zeopp import zeopp + + +@pytest.mark.parametrize("interface", ["legacy", "unified"]) +def test_unequal_radii_and_engine_output_names( + interface, sample_cif, tmp_path, monkeypatch +): + commands = [] + data = Path(__file__).parent / "data" / "zeopp" + + def engine(directory, stem, arguments): + commands.append(arguments) + for analysis in ("sa", "vol", "psd", "chan"): + suffix = "psd_histo" if analysis == "psd" else analysis + shutil.copyfile( + data / f"test_structure.{analysis}", + directory / f"{stem}.{suffix}", + ) + + if interface == "legacy": + monkeypatch.setattr(zeopp, "_find_network_binary", lambda _: "network") + + def launch(command, **kwargs): + cif = Path(command[-1]) + engine(cif.parent, cif.stem, command[1:]) + return SimpleNamespace(returncode=0, stderr="") + + monkeypatch.setattr(zeopp.subprocess, "run", launch) + result = zeopp.run_zeopp( + sample_cif, + analyses=["sa", "vol", "psd", "chan"], + probe_radius=1.2, + chan_radius=1.8, + num_samples=1234, + output_dir=str(tmp_path / "legacy"), + ) + assert result["success"] + results = result["results"] + else: + + def launch(root, execution, name, arguments): + engine(root / "work", "structure", arguments) + + monkeypatch.setattr(adapters, "external_command", launch) + result = run( + PoreRequest( + structure=StructureRef(path=sample_cif), + analyses=["sa", "vol", "psd", "chan"], + probe_radius=1.2, + channel_radius=1.8, + num_samples=1234, + ), + output_dir=tmp_path / "unified", + ) + assert result.accepted, result.failure + results = result.payload.results + + for flag in ("-sa", "-vol", "-psd"): + index = commands[0].index(flag) + assert commands[0][index + 1 : index + 4] == ["1.8", "1.2", "1234"] + index = commands[0].index("-chan") + assert commands[0][index + 1] == "1.2" + assert results["chan"]["dimensionalities"] == [3, 3] + assert results["psd"]["counts"][3] == 25 + + +@pytest.mark.parametrize("suffix", ["psd_histo", "psd"]) +def test_histogram_discovery_and_explicit_file(suffix, tmp_path): + path = tmp_path / f"structure.{suffix}" + path.write_text( + "# diameter count cumulative derivative\n0.1 2 0.5 1\n0.2 2 1 1\n" + ) + for target in (path, tmp_path): + result = zeopp.get_output_data(str(target), analyses=["psd"]) + assert result["success"] + assert result["psd"]["counts"] == [2, 2] + + +def test_zero_channels_and_header_only_legacy_file(tmp_path): + path = tmp_path / "structure.chan" + for count, dimensions in ((0, []), (1, [2])): + path.write_text( + f"structure {count} channels identified of dimensionality " + + " ".join(map(str, dimensions)) + + "\n" + ) + assert zeopp.get_output_data(str(path))["chan"] == { + "num_channels": count, + "dimensionalities": dimensions, + } + + +@pytest.mark.parametrize( + "contents", + [ + "structure 2 channels identified of dimensionality 3\n", + "structure 1 channels identified of dimensionality 4\n", + "structure 1 channels identified of dimensionality 3\n" + "Channel 0 1 nan 1\n", + "structure 1 channels identified of dimensionality 3\nChannel 0 1\n", + "structure 2 channels identified of dimensionality 3 3\n" + "Channel 0 1 1 1\n", + ], +) +def test_malformed_channels_fail(contents, tmp_path): + path = tmp_path / "structure.chan" + path.write_text(contents) + with pytest.raises(ValueError): + zeopp.get_output_data(str(path)) + + +@pytest.mark.parametrize( + "contents", + [ + "# header only\n", + "0.1\n", + "0.1 broken\n", + "0.1 nan\n", + "0.1 -1\n", + "0.1 1\n0.1 2\n", + "0.1 1\n0.2 2\n0.4 1\n", + ], +) +def test_malformed_histograms_fail(contents, tmp_path): + path = tmp_path / "structure.psd_histo" + path.write_text(contents) + with pytest.raises(ValueError): + zeopp.get_output_data(str(path)) + + +@pytest.mark.parametrize( + "analysis,contents", + [ + ("res", "structure nan 2 3\n"), + ("sa", "ASA_A^2: 1\n"), + ("vol", "AV_A^3: 1\n"), + ], +) +def test_incomplete_or_nonfinite_output_fails(analysis, contents, tmp_path): + path = tmp_path / f"structure.{analysis}" + path.write_text(contents) + with pytest.raises(ValueError): + zeopp.get_output_data(str(tmp_path), analyses=[analysis]) + + +def test_missing_requested_analysis_fails(sample_cif, tmp_path, monkeypatch): + monkeypatch.setattr(zeopp, "_find_network_binary", lambda _: "network") + monkeypatch.setattr( + zeopp.subprocess, + "run", + lambda *a, **k: SimpleNamespace(returncode=0, stderr=""), + ) + with pytest.raises(ValueError, match="Missing requested"): + zeopp.run_zeopp(sample_cif, output_dir=str(tmp_path / "run")) + with pytest.raises(ValueError, match="Missing requested"): + zeopp.get_output_data(str(tmp_path), analyses=["res"]) + + +def test_invalid_radius_relation_rejected_by_request(sample_cif): + with pytest.raises(ValueError, match="must not exceed"): + PoreRequest( + structure=StructureRef(path=sample_cif), + analyses=["sa"], + probe_radius=1.8, + channel_radius=1.2, + ) diff --git a/tests/test_zeopp.py b/tests/test_zeopp.py index 60ed90e..25f412f 100644 --- a/tests/test_zeopp.py +++ b/tests/test_zeopp.py @@ -19,6 +19,13 @@ ) +def _successful_res(cmd, **kwargs): + cif = Path(cmd[-1]) + source = Path(__file__).parent / "data/zeopp/test_structure.res" + shutil.copyfile(source, cif.with_suffix(".res")) + return MagicMock(returncode=0, stderr="") + + @pytest.fixture def zeopp_data_dir(test_data_dir): """Path to the zeopp test data directory.""" @@ -212,7 +219,7 @@ def test_run_includes_ha_flag(self, mock_run, mock_find, sample_cif, tmp_path): """Should include -ha flag by default.""" mock_find.return_value = "/usr/bin/network" - mock_run.return_value = MagicMock(returncode=0, stderr="") + mock_run.side_effect = _successful_res run_zeopp(sample_cif, analyses=["res"], output_dir=str(tmp_path / "out")) @@ -226,7 +233,7 @@ def test_run_uses_bundled_radii_by_default( ): """Should use bundled UFF.rad when no radii file given.""" mock_find.return_value = "/usr/bin/network" - mock_run.return_value = MagicMock(returncode=0, stderr="") + mock_run.side_effect = _successful_res outdir = tmp_path / "out" run_zeopp(sample_cif, analyses=["res"], @@ -240,7 +247,7 @@ def test_run_uses_bundled_radii_by_default( def test_run_no_ha_flag(self, mock_run, mock_find, sample_cif, tmp_path): """Should omit -ha flag when ha=False.""" mock_find.return_value = "/usr/bin/network" - mock_run.return_value = MagicMock(returncode=0, stderr="") + mock_run.side_effect = _successful_res run_zeopp(sample_cif, analyses=["res"], ha=False, output_dir=str(tmp_path / "out")) @@ -253,7 +260,7 @@ def test_run_with_radii_file(self, mock_run, mock_find, sample_cif, tmp_path): """Should include -r flag and copy radii file to workdir.""" mock_find.return_value = "/usr/bin/network" - mock_run.return_value = MagicMock(returncode=0, stderr="") + mock_run.side_effect = _successful_res # Create a fake radii file rad_file = tmp_path / "UFF.rad"