diff --git a/docs/cli-reference.md b/docs/cli-reference.md index 5fdf821..794e86b 100644 --- a/docs/cli-reference.md +++ b/docs/cli-reference.md @@ -127,8 +127,10 @@ compatibility. Use the three explicit selectors. | `--validation-engine` | `comunica`, `qlever`, `hdt`, `cottas`, `all`, or a comma-separated list | `comunica` | SPARQL backend(s); a scale and performance decision, never a semantic one. `hdt`/`cottas` query the compressed artifact in place. Several engines answer the whole query set, are cross-checked against each other, and are timed in `benchmark.csv` | | `--filter-oracle` | `auto`, `bcftools`, `cyvcf2` | `auto` | FILTER-field oracle | | `--validation-queries` | query ids, or `core` / `preflight` / `all` | all | Run only these queries. A subset reports `TIMING_ONLY` rather than a validation verdict, because the PASS decision needs the whole set; each selected query is still compared against the VCF oracle, so a disagreement still fails the run. Use it to measure retrieval cost for one question without paying for the rest | -| `--shacl-shapes` | path | off | Independent structural layer via `pyshacl`; in-memory, so not for cohort scale | -| `--shacl-max-triples` | N | 50,000,000 | Skip the shape layer, and record the skip, when the decoded graph exceeds N triples; `0` disables the gate. The gate reads the decoded graph, not the artifact it arrived in | +| `--shacl-shapes` | path | bundled `core` profile | Independent structural layer via `pyshacl`; replaces the bundled shapes | +| `--shacl-batch-triples` | N | 500,000 | Validate node-level shapes (the `core` profile) a batch of records of about N triples at a time, so memory follows the batch, not the graph; `0` validates whole. Shapes with SPARQL constraints are always validated whole | +| `--shacl-workers` | N | up to 4 | Shape batches validated in parallel; peak memory is about N batches | +| `--shacl-max-triples` | N | 10,000,000 | Skip shapes validated whole, and record the skip, when the decoded graph exceeds N triples; `0` disables the gate. The gate reads the decoded graph, not the artifact it arrived in | | `--node-heap-mb` | MB | Node's own | V8 old-space ceiling for the Comunica-backed engines (`comunica`, `hdt`, `cottas`); Node does not size its heap from the machine | | `--strict-conformance` | — | off | Promote a missing-token conformance anomaly from report to failure | | `--validation-query-timeout` | seconds | 3600 | Per-query timeout, every engine | diff --git a/docs/validation-methodology.md b/docs/validation-methodology.md index d6b5e59..20a483e 100644 --- a/docs/validation-methodology.md +++ b/docs/validation-methodology.md @@ -98,8 +98,9 @@ a plain literal instead of `xsd:decimal`/`vcfc:Null`; `##fileDate` untyped instead of `xsd:date`) and one of which is a contradiction inside the vocabulary itself, recorded in [`vcf-coverage.md`](vcf-coverage.md). -SHACL is opt-in because `pyshacl` loads the graph into memory and does not -scale to a cohort-sized aggregate. +`pyshacl` loads its data graph into memory, so node-level shapes are validated +a batch of records at a time (see [`validation.md`](validation.md#shacl)); +shapes with SPARQL constraints compare records and are validated whole. ### Identity digests, and why they are histograms @@ -113,6 +114,9 @@ own IRI**, then bucket on the first byte of the hash. Two properties matter: its order. Bucketing is order-independent by construction, and keeps the result at most 256 rows for a graph of any size. A mismatch is localized by re-querying only the differing buckets. +- **A digest compares values, not spellings.** QUAL is an `xsd:decimal`, which + QLever returns canonically (`30.1`) and other engines lexically (`30.10`), so + `q11` and its oracle both drop trailing fractional zeros before hashing. Fields are separated by U+001F, which cannot occur in a VCF field, so no shift of a field boundary can forge a match. diff --git a/docs/validation-migration-notes.md b/docs/validation-migration-notes.md index 21453f6..8a6745c 100644 --- a/docs/validation-migration-notes.md +++ b/docs/validation-migration-notes.md @@ -182,10 +182,10 @@ vocabulary repository ships one example VCF per version under if that coverage is wanted. **The SV carriers have no mutation.** They are represented and conformant, and -the census counts the families the fixture exercises, but the fixture contains -no structural variants, so `emitted_record_counters` is untested against -breakends, tandem repeats and reference blocks. Adding one SV record to the -fixture would exercise them through the existing machinery. +the census counts them -- `test_validation_real_files_unit.py` checks the +oracle's counts against the emitters' output for events, confidence intervals, +reference blocks and phase sets -- but no mutation targets them, and tandem +repeats are not counted yet. ## How conformance was verified diff --git a/docs/validation.md b/docs/validation.md index 96de7f2..8203f2e 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -575,8 +575,11 @@ The common record-level queries are used for both graph shapes: Q9 and Q10 are the completeness check: comparing the graph's inventory against what the VCF implies catches a predicate that is missing, one with the wrong -cardinality, and one that should not be there at all. Q11-Q13 close the -permutation gap - see +cardinality, and one that should not be there at all. The inventory covers phase +sets, the SV carriers (events, confidence intervals, IMPRECISE/NOVEL, SVLEN, +SVCLAIM) and gVCF reference blocks; tandem repeats, local alleles and base +modifications are not counted yet, so a file using them fails Q9/Q10. Q11-Q13 +close the permutation gap - see [`validation-methodology.md`](validation-methodology.md#identity-digests-and-why-they-are-histograms). Q9-Q13 assume the shipped RML mapping's predicate inventory and IRI templates. @@ -645,11 +648,21 @@ and a conforming graph is reported as violations. In a vocabulary checkout the bundle sits one level up from the shapes, in `ontology/`, and is found automatically; `--shacl-ontology PATH` names it explicitly anywhere else. -It is **on by default for sources at or below 512 MiB**, using shapes vendored -with the package -- no vocabulary checkout needed. `--no-shacl` turns it off; -`--shacl-shapes` overrides both the bundled shapes and the size gate. Above the -gate it is skipped, because `pyshacl` loads the whole graph into memory and a -cohort-scale aggregate would not fit. +It is **on by default**, using shapes vendored with the package -- no +vocabulary checkout needed. `--no-shacl` turns it off; `--shacl-shapes` +overrides the bundled shapes. + +`pyshacl` loads its data graph into memory, so node-level shapes -- the default +`core` profile -- are validated **a batch of records at a time**: each batch is +about `--shacl-batch-triples` triples (default 500,000) of whole records plus +the file-level triples they point at (header, sample set, definitions). Every +constraint in such a profile judges one node from its neighbourhood, so the +verdict is the whole-graph verdict; batches run in `--shacl-workers` processes +(default up to 4), and memory follows the batch, not the graph: about 1.3 GB +per default batch, so about 5 GB with four workers. Shapes with +SPARQL constraints compare records with each other, so they are validated whole +and stay under `--shacl-max-triples` (default 10,000,000), above which they are +skipped and the skip recorded. ### Two profiles, and which one catches what @@ -671,8 +684,9 @@ detects none of them. `full` is not the default because the cost is real and measured: its two extra profiles use `sh:sparql` constraints that self-join the graph, so they grow far -faster than the data. It is gated to sources at or below 16 MiB, where closing -those four classes is worth two minutes; the core profile's gate is 512 MiB. +faster than the data, and they compare records, so they are validated whole. +It is gated to sources at or below 16 MiB, where closing those four classes is +worth two minutes; the core profile has no size gate. ```bash # Close the four classes the query suite misses, on a fixture-sized input diff --git a/src/validation/queries/common/q11_record_digest.rq b/src/validation/queries/common/q11_record_digest.rq index d59e903..59c727a 100644 --- a/src/validation/queries/common/q11_record_digest.rq +++ b/src/validation/queries/common/q11_record_digest.rq @@ -15,6 +15,12 @@ PREFIX vcfc: # Fields are separated by U+001F (unit separator), written as a SPARQL UCHAR # escape. It cannot appear in a VCF field, so no combination of values can be # made to collide by shifting a field boundary. +# +# QUAL is an xsd:decimal, and engines disagree on how STR() spells one: QLever +# returns the canonical value ("30.1"), others the lexical form ("30.10"). So +# trailing zeros after the decimal point are dropped here, and the oracle drops +# them the same way (digest_qual), making the digest compare values rather than +# spellings. Nothing else changes: "100" stays "100" and "." stays ".". SELECT ?bucket (COUNT(*) AS ?recordCount) WHERE { ?record a vcfc:VCFRecord ; @@ -27,6 +33,9 @@ WHERE { ?call vcfc:qual ?qual ; vcfc:filter ?filter ; vcfc:infoRaw ?info . + BIND(STR(?qual) AS ?qualText) + BIND(IF(REGEX(?qualText, "\\.[0-9]*0$"), REPLACE(?qualText, "\\.?0+$", ""), ?qualText) + AS ?qualKey) BIND(SUBSTR(SHA256(CONCAT( STR(?record), "\u001F", STR(?chrom), "\u001F", @@ -34,7 +43,7 @@ WHERE { STR(?recordIdLiteral), "\u001F", STR(?ref), "\u001F", STR(?alt), "\u001F", - STR(?qual), "\u001F", + ?qualKey, "\u001F", STR(?filter), "\u001F", STR(?info) )), 1, 2) AS ?bucket) diff --git a/src/validation/validation_runner.py b/src/validation/validation_runner.py index 73d1005..3b87fc8 100644 --- a/src/validation/validation_runner.py +++ b/src/validation/validation_runner.py @@ -28,6 +28,7 @@ import gzip import hashlib import json +import multiprocessing import os import platform import re @@ -434,6 +435,11 @@ def classify_genotype(alleles: tuple[int | None, ...] | None, *, has_gt: bool) - RDF_TYPE = vocab.RDF_TYPE_URI +def _census_iri(name: str) -> str: + """A counter key as an IRI: a vcfc local name, or already a full IRI (FALDO).""" + return name if "://" in name else f"{VCFC}{name}" + + def _nonzero(counts: dict[str, int]) -> dict[str, int]: return {key: value for key, value in counts.items() if value} @@ -453,6 +459,18 @@ def _nonzero(counts: dict[str, int]) -> dict[str, int]: DIGEST_BUCKET_CHARS = 2 +def digest_qual(value: str) -> str: + """QUAL as q11_record_digest.rq hashes it: no trailing fractional zeros. + + The graph types QUAL as xsd:decimal, and engines spell a decimal differently + -- QLever canonically ("30.1"), others lexically ("30.10") -- so both sides + drop the zeros and the digest compares values. Integers and "." are kept. + """ + if re.search(r"\.[0-9]*0$", value): + return re.sub(r"\.?0+$", "", value) + return value + + def record_digest_bucket(fields: list[str]) -> str: """Bucket one record exactly as q11_record_digest.rq does.""" joined = DIGEST_SEPARATOR.join(fields) @@ -606,6 +624,106 @@ def _meta_values(raw: str | None) -> list[str]: return [item.strip() for item in inner.split(",") if item.strip()] +#: FALDO terms are counted by full IRI: CIPOS and CIEND intervals are +#: faldo:InRangePosition resources, not vcfc ones. +FALDO_IN_RANGE_POSITION = f"{vocab.FALDO_NAMESPACE}InRangePosition" +FALDO_BEGIN = f"{vocab.FALDO_NAMESPACE}begin" +FALDO_END = f"{vocab.FALDO_NAMESPACE}end" + + +def _count_sv_layer( + classes: Counter, predicates: Counter, *, alleles: list, info_map: dict, + pos: str, version: Any, +) -> None: + """Count the SV carriers ``_emit_sv_layer`` attaches to a record's ALT alleles. + + Every condition is the emitter's: a carrier exists only when the record + supplies what its shape needs. A repeated EVENT is one resource in the graph + however many alleles name it, so events are counted as distinct triples. + Tandem repeats are not modelled yet (see KNOWN_UNMODELLED in the tests). + """ + alts = [allele for allele in alleles if allele.index >= 1] + for key, predicate in vocab.SV_FLAG_INFO_KEYS.items(): + if key in info_map: + predicates[predicate] += len(alts) + for key, predicate in vocab.SV_ALLELE_INFO_KEYS.items(): + if info_map.get(key) is not None: + predicates[predicate] += min(len(vocab.split_value_items(info_map[key])), len(alts)) + if info_map.get("SVCLAIM") is not None: + claims = vocab.split_value_items(info_map["SVCLAIM"])[: len(alts)] + predicates["svClaim"] += sum( + 1 for claim in claims if claim.strip().upper() in vocab.SV_CLAIMS + ) + + # gVCF reference blocks: END on a <*>/ allele. + end = info_map.get("END") + if end is not None and end.lstrip("-").isdigit(): + blocks = sum(1 for allele in alts if allele.kind == "UnspecifiedAllele") + classes["ReferenceBlock"] += blocks + predicates["endPosition"] += blocks + predicates["isReferenceBlockStart"] += blocks + if pos.isdigit(): + predicates["referenceBlockLength"] += blocks + + event, event_type = info_map.get("EVENT"), info_map.get("EVENTTYPE") + if event and event_type and version.event_types: + names = vocab.split_value_items(event) + types = vocab.split_value_items(event_type) + if not version.events_per_alt: + names, types = names[:1], types[:1] + events, typed, links = set(), set(), set() + for offset, name in enumerate(names): + code = types[offset] if offset < len(types) else None + individual = vocab.EVENT_TYPES.get(code.strip().upper()) if code else None + if vocab.is_missing(name) or individual is None: + continue + if version.events_per_alt and offset >= len(alts): + continue + events.add(name) + typed.add((name, individual)) + links.add((offset + 1 if version.events_per_alt else 0, name)) + classes["VariantEvent"] += len(events) + predicates["eventType"] += len(typed) + predicates["inEvent"] += len(links) + + for key in (*vocab.FALDO_INTERVAL_INFO_KEYS, *vocab.CONFIDENCE_INTERVAL_INFO_KEYS): + if info_map.get(key) is None or version.tuple_arity(key) is None: + continue + items = vocab.split_value_items(info_map[key]) + for offset in range(len(alts)): + pair = offset if version.tuples_per_alt else 0 + if pair * 2 + 1 >= len(items): + break + if key in vocab.FALDO_INTERVAL_INFO_KEYS: + predicates[vocab.FALDO_INTERVAL_INFO_KEYS[key]] += 1 + classes[FALDO_IN_RANGE_POSITION] += 1 + predicates[FALDO_BEGIN] += 1 + predicates[FALDO_END] += 1 + else: + predicates[vocab.CONFIDENCE_INTERVAL_INFO_KEYS[key]] += 1 + classes["ConfidenceInterval"] += 1 + predicates["ciLower"] += 1 + predicates["ciUpper"] += 1 + + +def _count_phase_set(counts: Counter, fields: dict[str, str]) -> None: + """Count one sample call's phase set exactly as ``_emit_phase_set`` emits it.""" + identifier, name = fields.get("PS"), fields.get("PSL") + if vocab.is_missing(identifier) and vocab.is_missing(name): + return + counts["PhaseSet"] += 1 + counts["inPhaseSet"] += 1 + if not vocab.is_missing(identifier): + counts["phaseSetId"] += 1 + if not vocab.is_missing(name): + counts["phaseSetName"] += 1 + ordinal = fields.get("PSO") + if not vocab.is_missing(ordinal) and ordinal.lstrip("-").isdigit(): + counts["phaseSetOrdinal"] += 1 + if not vocab.is_missing(fields.get("PSQ")): + counts["phaseSetQuality"] += 1 + + def emitted_record_counters( rows: list[list[str]], samples: list[str], @@ -643,6 +761,7 @@ def emitted_record_counters( format_items = format_item_alleles = format_tuple_items = 0 format_item_predicates: Counter[str] = Counter() genotypes = genotype_calls = called_alleles = 0 + phase_set_counts: Counter[str] = Counter() for row in rows: chrom = row[0] if len(row) > 0 else "" @@ -705,6 +824,11 @@ def emitted_record_counters( predicates["forGenotypeIndex"] += 1 elif link.gt_allele_index is not None: predicates["forGTAlleleIndex"] += 1 + if entries and alt_count: + _count_sv_layer( + classes, predicates, alleles=alleles, info_map=dict(entries), + pos=row[1] if len(row) > 1 else "", version=version, + ) if samples: format_keys = (row[8].split(":") if len(row) > 8 and row[8] else []) @@ -746,6 +870,19 @@ def emitted_record_counters( format_item_predicates["forGenotypeIndex"] += 1 elif link.gt_allele_index is not None: format_item_predicates["forGTAlleleIndex"] += 1 + # A phase set is emitted per sample call from its non-empty PS/PSL/ + # PSO/PSQ cells, in the expanded profile only. + phase_keys = [ + (index, key) for index, key in enumerate(format_keys) + if key in vocab.PHASE_SET_FORMAT_KEYS + ] + if phase_keys: + for payload in payloads: + fields = payload.split(":") if payload else [] + _count_phase_set(phase_set_counts, { + key: fields[index] for index, key in phase_keys + if index < len(fields) and fields[index] + }) if "GT" in format_keys: gt_index = format_keys.index("GT") for payload in payloads: @@ -791,6 +928,10 @@ def emitted_record_counters( for name in ("hasAlleleCall", "callIndex", "isNoCall"): genotype_predicates[name] += genotype_calls genotype_predicates["calledAllele"] += called_alleles + # Phase sets travel with the expanded sample layer, as genotypes do. + if phase_set_counts: + genotype_classes["PhaseSet"] += phase_set_counts.pop("PhaseSet") + genotype_predicates.update(phase_set_counts) format_item_classes: Counter[str] = Counter() if format_items: @@ -1070,9 +1211,11 @@ def expected_census( # The value items and the SV carriers travel with the structured INFO # representation. The allele layer does not -- it is merged below. for class_name, count in parser["emittedRecordClasses"].items(): - classes[f"{VCFC}{class_name}"] = classes.get(f"{VCFC}{class_name}", 0) + count + iri = _census_iri(class_name) + classes[iri] = classes.get(iri, 0) + count for name, count in parser["emittedRecordPredicates"].items(): - predicates[f"{VCFC}{name}"] = predicates.get(f"{VCFC}{name}", 0) + count + iri = _census_iri(name) + predicates[iri] = predicates.get(iri, 0) + count # The allele layer is required by whoever joins to it: the structured INFO # value items (Number=A/R/G) or the expanded sample layer's calledAllele. @@ -1474,7 +1617,9 @@ def parse_vcf( f"file://{source_component}#record/" f"{rml_uri_component(str(total_records))}" ) - digest_buckets[record_digest_bucket([record_iri, *columns[:8]])] += 1 + digest_buckets[record_digest_bucket( + [record_iri, *columns[:5], digest_qual(columns[5]), *columns[6:8]] + )] += 1 # Structured INFO counts; `typed_value_kind` states the same rule # `_typed_field_object` applies in vcf_rdfizer.py. @@ -2024,18 +2169,19 @@ def merge_shapes_graph(shapes: list[Path]): return graph -#: pyshacl loads the whole graph into memory, so the shape layer needs a size -#: gate. This one is in triples rather than bytes, because what pyshacl pays -#: for is the graph, not the packaging -- see the comment at its call site for -#: the failure that motivated it. +#: pyshacl loads its whole data graph into memory, so shapes validated in one +#: piece -- those with SPARQL constraints, or with --shacl-batch-triples 0 -- +#: need a size gate; batched shapes do not. It is in triples rather than bytes, +#: because what pyshacl pays for is the graph, not the packaging -- see the +#: comment at its call site for the failure that motivated it. #: -#: 50M is chosen to sit above the largest graph the published campaign actually -#: validated with shapes (17.1M triples, the 100,000-record HG005 slice) and -#: below the one that exhausted 31 GB (170.9M). It is deliberately not derived -#: from measured bytes-per-triple: rdflib's footprint depends on term sharing -#: and IRI length, so a constant here is a conservative guard rather than a -#: prediction, and it is overridable for a machine that can afford more. -DEFAULT_SHACL_MAX_TRIPLES = 50_000_000 +#: Measured with the core shapes and RDFS inference on a real genome's graph, +#: pyshacl peaks at about 2.4 GB per million triples (1.26 GB at 0.5M, 4.72 GB +#: at 2.0M), so 10M is about 24 GB: it fits a 31 GB host, which 58.2M and +#: 170.9M did not. The published campaign's largest shape-validated graph was +#: 0.96M. rdflib's footprint depends on term sharing and IRI length, so this is +#: a conservative guard rather than a prediction; raise it on a larger machine. +DEFAULT_SHACL_MAX_TRIPLES = 10_000_000 #: Node's V8 heap ceiling for the Comunica-backed endpoints (comunica, hdt, #: cottas -- all three go through ComunicaHttpEndpointMixin). @@ -2082,11 +2228,125 @@ def shacl_exceeds_limit(triple_count: int | None, limit: int | None) -> bool: return triple_count > limit +#: Shapes are validated a batch of whole records at a time, so pyshacl's memory +#: is set by the batch rather than the graph. A batch holds its records' triples +#: and the file-level ones they point at (header, sample set, definitions), which +#: is everything a node-level constraint can see, so the verdict is the +#: whole-graph verdict. 0 validates the graph in one piece, as before. +#: +#: Time is linear in the batch (155 s at 0.5M triples, 607 s at 2.0M), so small +#: batches cost nothing in throughput; 0.5M peaks at about 1.3 GB. +DEFAULT_SHACL_BATCH_TRIPLES = 500_000 +#: Batches are independent and run in parallel; peak memory is about this many +#: batches at once, so about 5 GB at the defaults. +DEFAULT_SHACL_WORKERS = min(4, os.cpu_count() or 1) + +#: The record, its call and its sample calls -- and everything minted beneath +#: them -- carry the record's row in their IRI. Anything else is file-level. +RECORD_SCOPED_IRI = re.compile(r"#(?:record|call|sample)/([0-9]+)") +#: Lines held before write_shacl_batches appends them to their batch files. +SHACL_BATCH_FLUSH_LINES = 200_000 + + +def shapes_are_node_local(shapes: list[Path]) -> bool: + """Whether every constraint judges a node from its own neighbourhood. + + A SPARQL-based constraint can compare records with each other -- unique + record indices, nondecreasing POS -- which no batch can decide, so shapes + using one are validated in one piece and stay under the size gate. + """ + from rdflib.namespace import SH + + try: + graph = merge_shapes_graph([shapes] if isinstance(shapes, Path) else list(shapes)) + except Exception: # noqa: BLE001 - unparsable: validated whole, which reports it + return False + return (None, SH.sparql, None) not in graph + + +def shacl_records_per_batch( + records: int | None, triples: int | None, batch_triples: int, node_local: bool, +) -> int | None: + """Records per shape batch, or None to validate the graph in one piece.""" + if not (batch_triples and node_local and records and triples): + return None + if triples <= batch_triples: + return None + return max(1, batch_triples * records // triples) + + +def write_shacl_batches( + source: Path, workdir: Path, records_per_batch: int +) -> tuple[Path, list[Path]]: + """Split an N-Triples graph into record batches and one file-level context. + + A triple belongs to the record its subject is scoped to or, failing that, + its object's (vcfc:hasRecord points from the file to a record). Lines are + buffered and appended per batch, so no more than a few files are open. + """ + context_path = workdir / "context.nt" + buffers: dict[int, list[str]] = {} + buffered = 0 + + def flush() -> None: + for index, lines in buffers.items(): + with (workdir / f"batch-{index:06d}.nt").open("a", encoding="utf-8") as handle: + handle.writelines(lines) + buffers.clear() + + with source.open(encoding="utf-8") as lines, \ + context_path.open("w", encoding="utf-8") as context: + for line in lines: + subject, _, rest = line.partition(" ") + match = RECORD_SCOPED_IRI.search(subject) + if match is None: + obj = rest.partition(" ")[2] + if obj.startswith("<"): + match = RECORD_SCOPED_IRI.search(obj.partition(">")[0]) + if match is None: + context.write(line) + continue + buffers.setdefault((int(match.group(1)) - 1) // records_per_batch, []).append(line) + buffered += 1 + if buffered >= SHACL_BATCH_FLUSH_LINES: + flush() + buffered = 0 + flush() + return context_path, sorted(workdir.glob("batch-*.nt")) + + +def _validate_shacl_task(task: tuple) -> tuple[bool, str]: + """Run pyshacl over one data graph: a file-level context plus an optional batch.""" + from pyshacl import validate as pyshacl_validate + from rdflib import Graph + + data_path, batch_path, shapes, ontology = task + data = Graph() + data.parse(str(data_path), format="nt") + if batch_path is not None: + data.parse(str(batch_path), format="nt") + shacl_graph = str(shapes[0]) if len(shapes) == 1 else merge_shapes_graph(shapes) + conforms, _graph, text = pyshacl_validate( + data, + shacl_graph=shacl_graph, + ont_graph=str(ontology) if ontology is not None else None, + **({"shacl_graph_format": "turtle"} if isinstance(shacl_graph, str) else {}), + ont_graph_format="turtle" if ontology is not None else None, + inference="rdfs" if ontology is not None else "none", + advanced=True, + ) + return bool(conforms), text + + def validate_shacl( source: Path, shapes: Path | list[Path], results_dir: Path, ontology: Path | None = None, + *, + records_per_batch: int | None = None, + workers: int = DEFAULT_SHACL_WORKERS, + scratch_dir: Path | None = None, ) -> dict[str, Any]: """Validate the graph against SHACL shapes, if pyshacl is available. @@ -2101,14 +2361,16 @@ def validate_shacl( violations. Supplying it with RDFS inference is the configuration the vocabulary's own tests/validate_shacl.py uses. - pyshacl loads the graph into memory, so this is opt-in and unsuitable for a - cohort-scale aggregate. It is reported as EXECUTION_FAILED rather than a - conformance failure when the tool is missing, so an absent optional - dependency can never look like a bad graph. + pyshacl loads its data graph into memory, so with ``records_per_batch`` the + graph is validated a batch of records at a time (see write_shacl_batches), + in up to ``workers`` processes; results repeated across batches, such as a + file-level node's, are reported once. It is reported as EXECUTION_FAILED + rather than a conformance failure when the tool is missing, so an absent + optional dependency can never look like a bad graph. """ report_path = results_dir / "shacl.json" try: - from pyshacl import validate as pyshacl_validate + import pyshacl # noqa: F401 - presence check; the tasks import it except ImportError as error: result = { "status": "EXECUTION_FAILED", @@ -2120,22 +2382,25 @@ def validate_shacl( return result started = time.monotonic() + shapes_list = [shapes] if isinstance(shapes, Path) else list(shapes) try: - shapes_list = [shapes] if isinstance(shapes, Path) else list(shapes) - shacl_graph = ( - str(shapes_list[0]) if len(shapes_list) == 1 - else merge_shapes_graph(shapes_list) - ) - conforms, _graph, text = pyshacl_validate( - str(source), - shacl_graph=shacl_graph, - ont_graph=str(ontology) if ontology is not None else None, - data_graph_format="nt", - **({"shacl_graph_format": "turtle"} if isinstance(shacl_graph, str) else {}), - ont_graph_format="turtle" if ontology is not None else None, - inference="rdfs" if ontology is not None else "none", - advanced=True, - ) + with tempfile.TemporaryDirectory( + dir=scratch_dir or results_dir, prefix="shacl-batches-" + ) as work: + if records_per_batch: + context, batches = write_shacl_batches(source, Path(work), records_per_batch) + tasks = [(context, batch, shapes_list, ontology) for batch in batches] + tasks = tasks or [(context, None, shapes_list, ontology)] + else: + tasks = [(source, None, shapes_list, ontology)] + workers = max(1, min(workers, len(tasks))) + if workers == 1: + outcomes = [_validate_shacl_task(task) for task in tasks] + else: + # One batch per worker process: rdflib does not hand freed + # memory back, so a reused worker would hold its largest batch. + with multiprocessing.get_context("fork").Pool(workers, maxtasksperchild=1) as pool: + outcomes = pool.map(_validate_shacl_task, tasks) except Exception as error: # noqa: BLE001 - reported, never fatal here result = { "status": "EXECUTION_FAILED", @@ -2146,7 +2411,16 @@ def validate_shacl( write_json(report_path, result) return result - results = parse_shacl_results(text) + conforms = all(ok for ok, _text in outcomes) + text = "\n\n".join(text for _ok, text in outcomes) + # Each batch's report is parsed on its own, so no block runs into the next + # report's header, then results repeated across batches are kept once. + results, seen = [], set() + for _ok, batch_text in outcomes: + for entry in parse_shacl_results(batch_text): + if entry["text"] not in seen: + seen.add(entry["text"]) + results.append(entry) # Unknown counts as blocking: a result this parser could not classify is a # parser bug, and the safe reading of a parser bug is not "conformant". violations = [ @@ -2180,6 +2454,9 @@ def validate_shacl( "advisorySample": [entry["text"] for entry in advisories][:SHACL_SAMPLE_LIMIT], "report": str(log_path), "wallSeconds": time.monotonic() - started, + "batches": len(tasks), + "recordsPerBatch": records_per_batch, + "workers": workers, "sampleLimitedTo": SHACL_SAMPLE_LIMIT, "sample": violations[:SHACL_SAMPLE_LIMIT], } @@ -4043,8 +4320,17 @@ def run_validation(args: argparse.Namespace) -> int: shacl_limit = getattr( args, "shacl_max_triples", DEFAULT_SHACL_MAX_TRIPLES ) - if args.shacl_shapes is not None and shacl_exceeds_limit( - decoded_triples, shacl_limit + records_per_batch = None + if args.shacl_shapes is not None: + records_per_batch = shacl_records_per_batch( + parser.get("totalRecords"), decoded_triples, + getattr(args, "shacl_batch_triples", DEFAULT_SHACL_BATCH_TRIPLES), + shapes_are_node_local(args.shacl_shapes), + ) + # A batched validation's memory is set by its batch, not the graph, + # so only a graph validated in one piece is gated. + if args.shacl_shapes is not None and records_per_batch is None and ( + shacl_exceeds_limit(decoded_triples, shacl_limit) ): # An unknown count is refused too (see shacl_exceeds_limit), so # the reason must not assume there is a number to format. @@ -4073,7 +4359,10 @@ def run_validation(args: argparse.Namespace) -> int: if args.shacl_shapes is not None and shacl_skipped is None: progress.emit("progress", completed=0, detail="validating SHACL shapes") shacl_result = validate_shacl( - decoded, args.shacl_shapes, results_dir, args.shacl_ontology + decoded, args.shacl_shapes, results_dir, args.shacl_ontology, + records_per_batch=records_per_batch, + workers=getattr(args, "shacl_workers", DEFAULT_SHACL_WORKERS), + scratch_dir=getattr(args, "scratch_dir", None), ) if shacl_skipped is not None: shacl_result = shacl_skipped @@ -4479,15 +4768,33 @@ def build_arg_parser() -> argparse.ArgumentParser: "keeps Node's default" ), ) + parser.add_argument( + "--shacl-batch-triples", + type=int, + default=DEFAULT_SHACL_BATCH_TRIPLES, + help=( + "Validate node-level shapes a batch of records at a time, about " + "this many triples per batch, so memory follows the batch rather " + "than the graph. Shapes with SPARQL constraints compare records and " + f"are validated whole (0 always validates whole; default: " + f"{DEFAULT_SHACL_BATCH_TRIPLES:,})" + ), + ) + parser.add_argument( + "--shacl-workers", + type=int, + default=DEFAULT_SHACL_WORKERS, + help=f"Shape batches validated in parallel (default: {DEFAULT_SHACL_WORKERS})", + ) parser.add_argument( "--shacl-max-triples", type=int, default=DEFAULT_SHACL_MAX_TRIPLES, help=( - "Skip the shape layer when the decoded graph holds more triples " - "than this, recording the skip and its reason (0 disables the " - "gate). pyshacl is in-memory, and its cost tracks the graph rather " - f"than the artifact it arrived in (default: {DEFAULT_SHACL_MAX_TRIPLES:,})" + "Skip shapes validated whole (see --shacl-batch-triples) when the " + "decoded graph holds more triples than this, recording the skip and " + "its reason (0 disables the gate). pyshacl is in-memory, and its cost " + f"tracks the graph, not the artifact (default: {DEFAULT_SHACL_MAX_TRIPLES:,})" ), ) parser.add_argument("--filter-oracle", choices=("auto", "bcftools", "cyvcf2"), default="auto") diff --git a/test/test_shacl_default_unit.py b/test/test_shacl_default_unit.py index 235a14d..03cffbc 100644 --- a/test/test_shacl_default_unit.py +++ b/test/test_shacl_default_unit.py @@ -77,12 +77,8 @@ def test_the_full_profile_resolves_all_three_files(self): ) self.assertIsNotNone(ontology) - def test_the_full_profile_has_a_much_smaller_size_gate(self): - """It is quadratic-ish in graph size; the core gate would be unusable.""" - self.assertLess( - vcf_rdfizer.FULL_SHACL_MAX_SOURCE_BYTES, - vcf_rdfizer.DEFAULT_SHACL_MAX_SOURCE_BYTES, - ) + def test_only_the_full_profile_is_size_gated(self): + """It is validated whole and quadratic-ish; core is validated in batches.""" big = vcf_rdfizer.FULL_SHACL_MAX_SOURCE_BYTES + 1 self.assertFalse(vcf_rdfizer.shacl_default_applies(big, "full")) self.assertTrue(vcf_rdfizer.shacl_default_applies(big, "core")) @@ -162,31 +158,20 @@ def test_the_vendored_copies_match_the_vocabulary_checkout(self): class SizeGateTests(VerboseTestCase): - def test_a_small_source_gets_shapes_by_default(self): - self.assertTrue(vcf_rdfizer.shacl_default_applies(10 * 1024 * 1024)) - - def test_a_source_at_the_limit_still_gets_shapes(self): - self.assertTrue( - vcf_rdfizer.shacl_default_applies( - vcf_rdfizer.DEFAULT_SHACL_MAX_SOURCE_BYTES - ) - ) - - def test_a_cohort_scale_source_does_not(self): - """pyshacl is in-memory; trying anyway turns a safety net into an OOM.""" - self.assertFalse( - vcf_rdfizer.shacl_default_applies( - vcf_rdfizer.DEFAULT_SHACL_MAX_SOURCE_BYTES + 1 - ) - ) - self.assertFalse(vcf_rdfizer.shacl_default_applies(200 * 1024 ** 3)) - - def test_an_unknown_size_is_treated_as_too_large(self): + def test_the_core_profile_applies_at_any_size(self): + """The validator batches it by record, so its memory does not follow the input.""" + for size in (10 * 1024 * 1024, 200 * 1024 ** 3, None): + with self.subTest(size=size): + self.assertTrue(vcf_rdfizer.shacl_default_applies(size)) + + def test_the_full_profile_applies_up_to_its_gate(self): + limit = vcf_rdfizer.FULL_SHACL_MAX_SOURCE_BYTES + self.assertTrue(vcf_rdfizer.shacl_default_applies(limit, "full")) + self.assertFalse(vcf_rdfizer.shacl_default_applies(limit + 1, "full")) + + def test_an_unknown_size_is_too_large_for_the_full_profile(self): """Skipping a check is recoverable; exhausting memory mid-run is not.""" - self.assertFalse(vcf_rdfizer.shacl_default_applies(None)) - - def test_the_gate_is_a_documented_constant_not_a_literal(self): - self.assertGreater(vcf_rdfizer.DEFAULT_SHACL_MAX_SOURCE_BYTES, 0) + self.assertFalse(vcf_rdfizer.shacl_default_applies(None, "full")) class CliTests(VerboseTestCase): diff --git a/test/test_shacl_gate_and_node_heap_unit.py b/test/test_shacl_gate_and_node_heap_unit.py index da7dc79..a9dfab1 100644 --- a/test/test_shacl_gate_and_node_heap_unit.py +++ b/test/test_shacl_gate_and_node_heap_unit.py @@ -86,9 +86,13 @@ def test_the_graph_that_oomed_is_refused(self): self.assertTrue(V.shacl_exceeds_limit(GRAPH_TRIPLES, V.DEFAULT_SHACL_MAX_TRIPLES)) def test_the_campaigns_largest_validated_graph_still_passes(self): - """17.1M triples is the 100,000-record HG005 slice the campaign ran - shapes on. The gate must not retroactively disable published behaviour.""" - self.assertFalse(V.shacl_exceeds_limit(17_098_746, V.DEFAULT_SHACL_MAX_TRIPLES)) + """0.96M triples is the largest graph the campaign ran shapes on. The + gate must not retroactively disable published behaviour.""" + self.assertFalse(V.shacl_exceeds_limit(958_919, V.DEFAULT_SHACL_MAX_TRIPLES)) + + def test_the_real_genome_that_exhausted_memory_is_refused(self): + """58.2M triples of NG131FQA1I hung a 31 GB host in one piece.""" + self.assertTrue(V.shacl_exceeds_limit(58_231_176, V.DEFAULT_SHACL_MAX_TRIPLES)) def test_an_unknown_count_is_treated_as_too_large(self): """Skipping is recoverable and recorded; exhausting memory is not.""" @@ -107,8 +111,8 @@ def test_the_boundary_is_inclusive(self): def test_the_default_sits_between_the_two_observations(self): """It must admit what worked and refuse what died, or it is arbitrary.""" - self.assertGreater(V.DEFAULT_SHACL_MAX_TRIPLES, 17_098_746) - self.assertLess(V.DEFAULT_SHACL_MAX_TRIPLES, GRAPH_TRIPLES) + self.assertGreater(V.DEFAULT_SHACL_MAX_TRIPLES, 958_919) + self.assertLess(V.DEFAULT_SHACL_MAX_TRIPLES, 58_231_176) class NodeHeapEnvTests(VerboseTestCase): @@ -271,6 +275,9 @@ def _runner_args(tmp: Path, **overrides) -> argparse.Namespace: qlever_server_arg=[], vcf=tmp / "s.vcf", filter_oracle="cyvcf2", dataset_id="sample", queries=None, shacl_max_triples=V.DEFAULT_SHACL_MAX_TRIPLES, + # Whole-graph validation, the path the size gate guards. The batched + # path is exempt from it; RunValidationShaclGateTests pins both. + shacl_batch_triples=0, shacl_workers=1, node_heap_mb=V.DEFAULT_NODE_HEAP_MB, ) fields.update(overrides) @@ -364,6 +371,26 @@ def test_a_graph_above_the_limit_never_reaches_pyshacl(self): self.assertRunDidNotError(run) run["validate_shacl"].assert_not_called() + def test_batched_shapes_are_not_gated(self): + """A batch bounds pyshacl's memory, so a large graph still gets shapes.""" + run = self._drive( + rdf_validation={"status": "PASS", "tripleCount": GRAPH_TRIPLES}, + shacl_batch_triples=V.DEFAULT_SHACL_BATCH_TRIPLES, shacl_workers=3) + self.assertRunDidNotError(run) + run["validate_shacl"].assert_called_once() + self.assertTrue(run["validate_shacl"].call_args.kwargs["records_per_batch"]) + self.assertEqual(run["validate_shacl"].call_args.kwargs["workers"], 3) + self.assertIsNone(run["shacl"], "nothing was skipped, so nothing is recorded") + + def test_sparql_shapes_are_still_gated_when_batching_is_on(self): + """A SPARQL constraint compares records, so those shapes stay whole.""" + with mock.patch.object(V, "shapes_are_node_local", return_value=False): + run = self._drive( + rdf_validation={"status": "PASS", "tripleCount": GRAPH_TRIPLES}, + shacl_batch_triples=V.DEFAULT_SHACL_BATCH_TRIPLES) + run["validate_shacl"].assert_not_called() + self.assertEqual(run["shacl"]["status"], "SKIPPED_TOO_LARGE") + def test_the_skip_is_recorded_with_its_reason(self): """A skip nobody can see is the silent failure this replaces.""" run = self._drive( @@ -398,7 +425,7 @@ def test_a_skip_is_not_a_shape_failure(self): def test_a_graph_within_the_limit_still_gets_its_shapes_checked(self): """The campaign's largest shape-validated graph, at the default gate.""" run = self._drive( - rdf_validation={"status": "PASS", "tripleCount": 17_098_746}) + rdf_validation={"status": "PASS", "tripleCount": 958_919}) self.assertRunDidNotError(run) run["validate_shacl"].assert_called_once() self.assertEqual(run["validate_shacl"].call_args[0][0], run["decoded"]) @@ -647,6 +674,24 @@ def test_a_bad_triple_limit_is_a_usage_error(self): self.assertIsNone(options, "a bad value must stop before the run") self.assertIn(message, err) + def test_the_batch_options_are_forwarded_as_integers(self): + _rc, options, _err = self._main( + ["--shacl-batch-triples", "0", "--shacl-workers", "2"]) + self.assertEqual(options["shacl_batch_triples"], 0) + self.assertEqual(options["shacl_workers"], 2) + + def test_bad_batch_options_are_usage_errors(self): + for argv, message in ( + (["--shacl-batch-triples", "-1"], + "--shacl-batch-triples must be zero or a positive integer"), + (["--shacl-workers", "0"], "--shacl-workers must be a positive integer"), + ): + with self.subTest(argv=argv): + rc, options, err = self._main(argv) + self.assertEqual(rc, 2) + self.assertIsNone(options) + self.assertIn(message, err) + def test_the_heap_is_forwarded_as_an_integer(self): _rc, options, _err = self._main(["--node-heap-mb", "16384"]) self.assertEqual(options["node_heap_mb"], 16384) @@ -712,9 +757,11 @@ def test_the_runner_accepts_what_the_wrapper_sends(self): """The contract across the boundary: the flag names and value types the wrapper emits are ones the runner parses.""" command = self._command_for( - {"shacl_max_triples": 0, "node_heap_mb": 16384}) + {"shacl_max_triples": 0, "node_heap_mb": 16384, + "shacl_batch_triples": 0, "shacl_workers": 2}) actions = V.build_arg_parser()._option_string_actions - for flag, expected in (("--shacl-max-triples", 0), ("--node-heap-mb", 16384)): + for flag, expected in (("--shacl-max-triples", 0), ("--node-heap-mb", 16384), + ("--shacl-batch-triples", 0), ("--shacl-workers", 2)): with self.subTest(flag=flag): self.assertIn(flag, actions) self.assertEqual(actions[flag].type(self._value(command, flag)), expected) diff --git a/test/test_validation_oracle_unit.py b/test/test_validation_oracle_unit.py index 581a48c..5473c6a 100644 --- a/test/test_validation_oracle_unit.py +++ b/test/test_validation_oracle_unit.py @@ -885,14 +885,6 @@ class EmittedTermCensusCoverageTests(VerboseTestCase): # VCF 4.5 local alleles (LA/LR/LG) "LocalAlleleSet", "LocalAlleleMembership", "hasLocalAlleleSet", "hasLocalAlleleMembership", "hasLocalAllele", "localAllele", "localIndex", - # phase sets (PS/PSL/PSO/PSQ) - "PhaseSet", "inPhaseSet", "phaseSetId", "phaseSetName", "phaseSetOrdinal", - "phaseSetQuality", - # structural variants and confidence intervals - "VariantEvent", "inEvent", "eventType", "svClaim", - "ConfidenceInterval", "ciLower", "ciUpper", "endPosition", - # gVCF reference blocks - "ReferenceBlock", "isReferenceBlockStart", "referenceBlockLength", # tandem repeats "TandemRepeatAllele", "RepeatSequence", "hasRepeatSequence", "repeatSequenceCount", "repeatSequenceIndex", diff --git a/test/test_validation_real_files_unit.py b/test/test_validation_real_files_unit.py new file mode 100644 index 0000000..7974c8f --- /dev/null +++ b/test/test_validation_real_files_unit.py @@ -0,0 +1,395 @@ +"""Validation on real files: three gaps a real genome exposed. + +Paired validation of 250,000 records of a real GATK genome (NG131FQA1I, 58.2M +triples) reported MISMATCH although every field of every record equals the VCF +text. Two of the causes were the validator's, and a third made the run itself +impossible: + +* **The census did not count phase sets** (30,910 from GATK's PS), nor the SV + and gVCF carriers -- reference blocks, events, confidence intervals, IMPRECISE, + SVLEN -- so any file using them failed with extra rows. +* **The record digest hashed STR(QUAL)**, which QLever returns canonically + ("30.1") and other engines lexically ("30.10"): 1,998 records hashed apart. +* **pyshacl loaded the whole graph**, and the 277 MB artifact exhausted a 31 GB + host. Node-level shapes are now validated a batch of records at a time. +""" + +import importlib.util +import os +import sys +import tempfile +import types +import unittest +from collections import Counter +from pathlib import Path +from unittest import mock + +import vcf_rdfizer +import vcf_rdfizer_vocab as vocab +from test.helpers import VerboseTestCase + +RUNNER_PATH = Path(__file__).resolve().parents[1] / "src" / "validation" / "validation_runner.py" +_spec = importlib.util.spec_from_file_location("validation_runner_real_files", RUNNER_PATH) +V = importlib.util.module_from_spec(_spec) +# Registered so the parallel shape workers can unpickle their task function. +sys.modules[_spec.name] = V +_spec.loader.exec_module(V) + +VCFC = vocab.VCFC_NAMESPACE +RDF_TYPE = vocab.RDF_TYPE_URI +FALDO = vocab.FALDO_NAMESPACE + +try: + import pyshacl # noqa: F401 + HAVE_PYSHACL = True +except ImportError: + HAVE_PYSHACL = False + + +def _emit(tmp: Path, file_format: str, rows: list[list[str]], headers: list[tuple[str, str]]): + """Run the expanded-sample and structured-INFO emitters the way full mode does.""" + records = tmp / "s.records.tsv" + records.write_text( + "SOURCE_FILE\tROW_ID\tCHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS1\n" + + "".join(f"s.vcf\t{i}\t" + "\t".join(row) + "\n" for i, row in enumerate(rows, 1)), + encoding="utf-8", + ) + lines = [("fileformat", file_format), *headers] + header_tsv = tmp / "s.header_lines.tsv" + header_tsv.write_text( + "SOURCE_FILE\tHEADER_INDEX\tHEADER_KEY\tHEADER_VALUE\tRAW_LINE\n" + + "".join(f"s.vcf\t{i}\t{key}\t{value}\tx\n" for i, (key, value) in enumerate(lines, 1)), + encoding="utf-8", + ) + rdf = tmp / "s.nt" + rdf.write_text("", encoding="utf-8") + # The wrapper resolves the version from the header and passes it on; the + # emitters fall back to the newest version otherwise. + version, _ = vocab.resolve_vcf_version(file_format) + vcf_rdfizer.append_expanded_sample_rdf( + records, rdf, header_tsv, progress_interval_records=0, version=version) + vcf_rdfizer.emit_record_detail( + "structured", records_tsv=records, header_lines_tsv=header_tsv, rdf_path=rdf, + sample_representation="expanded", version=version, + ) + return rdf + + +class SvGvcfPhaseSetCensusTests(VerboseTestCase): + """The oracle counts what the emitters write, term by term, for every family.""" + + HEADERS = [ + ("contig", ""), + ("INFO", ""), + ("INFO", ""), + ("INFO", ""), + ("FORMAT", ""), + ("FORMAT", ""), + ] + ROWS = [ + # Two SV alleles in one event: the event is one resource with two links. + ["chr1", "100", ".", "A", ",", "50", "PASS", + "SVLEN=-100,50;IMPRECISE;CIPOS=-5,5,-3,3;CIEND=-2,2,-1,1;CILEN=-10,10,-4,4;" + "SVCLAIM=D,J;EVENT=ev1,ev1;EVENTTYPE=DEL,DEL", "GT:PS", "0/1:100"], + # gVCF reference blocks, alone and beside an explicit ALT. + ["chr1", "200", ".", "G", "<*>", "0", ".", "END=250", "GT:PS:PSO", "0/0:.:."], + ["chr1", "400", ".", "T", "C,", "20", ".", "END=450", "GT", "0/1"], + # Phase sets: complete, name-only, and missing (no carrier). + ["chr1", "300", ".", "C", "T", "40", "PASS", "DP=5", + "GT:PS:PSL:PSO:PSQ", "0|1:300:blockA:2:30"], + ["chr1", "310", ".", "C", "G", "40", "PASS", "DP=5", "GT:PS:PSL", "0|1:.:blockB"], + ["chr1", "320", ".", "C", "A", "40", "PASS", "DP=5", "GT:PS", "0|1:."], + # SV INFO on a record with no ALT: there is no allele to carry it. + ["chr1", "500", ".", "A", ".", "10", "PASS", "SVLEN=5;IMPRECISE;CIPOS=-1,1", "GT", "0/0"], + ] + PREDICATES = { + *(VCFC + name for name in ( + "isImprecise", "isNovel", "svLength", "svClaim", "endPosition", + "referenceBlockLength", "isReferenceBlockStart", "inEvent", "eventType", + "posConfidenceInterval", "endConfidenceInterval", "lenConfidenceInterval", + "copyNumberConfidenceInterval", "ciLower", "ciUpper", "inPhaseSet", + "phaseSetId", "phaseSetName", "phaseSetOrdinal", "phaseSetQuality")), + FALDO + "begin", FALDO + "end", + } + CLASSES = { + *(VCFC + name for name in ( + "ReferenceBlock", "VariantEvent", "ConfidenceInterval", "PhaseSet")), + FALDO + "InRangePosition", + } + + def emitted(self, file_format): + """Per-term counts in the emitted graph, as a set of triples.""" + with tempfile.TemporaryDirectory() as td: + lines = set(_emit(Path(td), file_format, self.ROWS, self.HEADERS) + .read_text(encoding="utf-8").splitlines()) + predicates, classes = Counter(), Counter() + for line in lines: + parts = line.split(" ", 2) + if len(parts) < 3: + continue + predicate, obj = parts[1][1:-1], parts[2].rsplit(" .", 1)[0] + if predicate in self.PREDICATES: + predicates[predicate] += 1 + if predicate == RDF_TYPE and obj[1:-1] in self.CLASSES: + classes[obj[1:-1]] += 1 + return predicates, classes + + def counted(self, file_format): + """The same terms as the oracle counts them from the VCF columns.""" + version, _ = vocab.resolve_vcf_version(file_format) + out = V.emitted_record_counters( + [list(row) for row in self.ROWS], ["S1"], version=version, + contig_ids={"chr1"}, alt_declaration_ids=set(), + info_numbers={"SVLEN": "A", "END": "1", "CIPOS": "."}, + format_numbers={"GT": "1", "PS": "1"}, + ) + predicates, classes = Counter(), Counter() + for bucket, target, terms in ( + ("emittedRecordPredicates", predicates, self.PREDICATES), + ("emittedGenotypePredicates", predicates, self.PREDICATES), + ("emittedRecordClasses", classes, self.CLASSES), + ("emittedGenotypeClasses", classes, self.CLASSES), + ): + for name, count in out[bucket].items(): + if V._census_iri(name) in terms and count: + target[V._census_iri(name)] += count + return predicates, classes + + def test_oracle_and_emitter_agree_on_every_term(self): + for file_format in ("VCFv4.2", "VCFv4.5"): + with self.subTest(version=file_format): + emitted_predicates, emitted_classes = self.emitted(file_format) + counted_predicates, counted_classes = self.counted(file_format) + self.assertEqual(dict(counted_predicates), dict(emitted_predicates)) + self.assertEqual(dict(counted_classes), dict(emitted_classes)) + + def test_the_fixture_exercises_every_family(self): + """A family the fixture never emits would agree vacuously.""" + predicates, classes = self.emitted("VCFv4.5") + for term in ("PhaseSet", "ReferenceBlock", "VariantEvent", "ConfidenceInterval"): + self.assertIn(VCFC + term, classes, term) + self.assertIn(FALDO + "InRangePosition", classes) + for term in ("isImprecise", "svLength", "svClaim", "phaseSetName", "phaseSetOrdinal"): + self.assertIn(VCFC + term, predicates, term) + + def test_a_shared_event_is_one_resource(self): + _predicates, classes = self.counted("VCFv4.5") + self.assertEqual(classes[VCFC + "VariantEvent"], 1) + + def test_events_need_vcf_4_4(self): + """Before 4.4 there is no EVENTTYPE, so EVENT stays an ordinary value.""" + _predicates, classes = self.counted("VCFv4.2") + self.assertNotIn(VCFC + "VariantEvent", classes) + + def test_counter_keys_become_iris(self): + self.assertEqual(V._census_iri("inPhaseSet"), VCFC + "inPhaseSet") + self.assertEqual(V._census_iri(FALDO + "begin"), FALDO + "begin") + + +class QualDigestTests(VerboseTestCase): + """Q11 compares QUAL values, not the spelling an engine returns.""" + + def test_trailing_fractional_zeros_are_dropped(self): + for value, expected in ( + ("30.10", "30.1"), ("100.0", "100"), ("100.00", "100"), ("0.50", "0.5"), + ("100", "100"), ("10", "10"), (".", "."), ("30.1", "30.1"), ("1e10", "1e10"), + ): + with self.subTest(value=value): + self.assertEqual(V.digest_qual(value), expected) + + def test_the_query_lands_in_the_oracle_bucket_for_every_spelling(self): + """The real Q11 text, run by rdflib's engine, against the oracle's hash.""" + from rdflib import Graph + + query = (V.QUERY_ROOT / "common" / "q11_record_digest.rq").read_text(encoding="utf-8") + record, call = "file://s.vcf#record/1", "file://s.vcf#call/1" + for qual in ("30.10", "30.1", "100.0", "100", "7"): + with self.subTest(qual=qual): + graph = Graph() + graph.parse(data="\n".join([ + f"<{record}> <{RDF_TYPE}> <{VCFC}VCFRecord> .", + f'<{record}> <{VCFC}chrom> "1" .', + f'<{record}> <{VCFC}pos> "5"^^ .', + f'<{record}> <{VCFC}recordId> "." .', + f'<{record}> <{VCFC}ref> "A" .', + f'<{record}> <{VCFC}alt> "G" .', + f"<{record}> <{VCFC}hasCall> <{call}> .", + f'<{call}> <{VCFC}qual> "{qual}"^^ .', + f'<{call}> <{VCFC}filter> "PASS" .', + f'<{call}> <{VCFC}infoRaw> "." .', + ]), format="nt") + buckets = [str(row.bucket) for row in graph.query(query)] + expected = V.record_digest_bucket( + [record, "1", "5", ".", "A", "G", V.digest_qual(qual), "PASS", "."]) + self.assertEqual(buckets, [expected]) + + +class ShaclBatchTests(VerboseTestCase): + """Node-level shapes are validated in record batches with the same verdict.""" + + GRAPH = "\n".join([ + " .", + " \"S1\" .", + *(f" <{VCFC}hasRecord> ." for n in range(1, 6)), + *(f" ." for n in range(1, 6)), + *(f" ." for n in range(1, 6)), + ]) + "\n" + + def test_every_triple_goes_to_exactly_one_place(self): + with tempfile.TemporaryDirectory() as td: + tmp = Path(td) + source = tmp / "g.nt" + source.write_text(self.GRAPH, encoding="utf-8") + context, batches = V.write_shacl_batches(source, tmp, records_per_batch=2) + parts = [context.read_text().splitlines()] + [b.read_text().splitlines() for b in batches] + self.assertEqual(len(batches), 3) + self.assertEqual(sorted(sum(parts, [])), sorted(self.GRAPH.splitlines())) + self.assertEqual(len(parts[0]), 2, "only the file-level triples are context") + # hasRecord follows its object: records 1-2 land in the first batch. + self.assertTrue(any("hasRecord> " in line for line in parts[1])) + self.assertFalse(any("#record/3" in line or "#sample/3/" in line for line in parts[1])) + + def test_batching_is_chosen_only_where_it_is_exact_and_needed(self): + rpb = V.shacl_records_per_batch + self.assertEqual(rpb(1000, 10_000_000, 2_000_000, True), 200) + self.assertEqual(rpb(10, 10_000_000, 2_000_000, True), 2) + self.assertEqual(rpb(10, 100_000_000, 2_000_000, True), 1) + self.assertIsNone(rpb(1000, 1_000_000, 2_000_000, True), "fits one batch") + self.assertIsNone(rpb(1000, 10_000_000, 0, True), "0 disables batching") + self.assertIsNone(rpb(1000, 10_000_000, 2_000_000, False), "SPARQL shapes") + self.assertIsNone(rpb(None, 10_000_000, 2_000_000, True)) + self.assertIsNone(rpb(1000, None, 2_000_000, True)) + # The real genome at the default: the 117 batches of its end-to-end run. + self.assertEqual(rpb(250_000, 58_231_176, V.DEFAULT_SHACL_BATCH_TRIPLES, True), 2146) + + def test_buffered_lines_are_appended_not_overwritten(self): + """Flushing after every line must give the same batches as one flush.""" + def split(flush_lines): + with tempfile.TemporaryDirectory() as td, \ + mock.patch.object(V, "SHACL_BATCH_FLUSH_LINES", flush_lines): + tmp = Path(td) + source = tmp / "g.nt" + source.write_text(self.GRAPH, encoding="utf-8") + context, batches = V.write_shacl_batches(source, tmp, records_per_batch=2) + return context.read_text(), [b.read_text() for b in batches] + + self.assertEqual(split(1), split(V.SHACL_BATCH_FLUSH_LINES)) + + def test_only_sparql_free_shapes_are_node_local(self): + core, _ = vcf_rdfizer.resolve_default_shacl_shapes(Path(vcf_rdfizer.__file__).parent, "core") + full, _ = vcf_rdfizer.resolve_default_shacl_shapes(Path(vcf_rdfizer.__file__).parent, "full") + self.assertTrue(V.shapes_are_node_local(core)) + self.assertFalse(V.shapes_are_node_local(full)) + with tempfile.TemporaryDirectory() as td: + broken = Path(td) / "broken.ttl" + broken.write_text("this is not turtle", encoding="utf-8") + self.assertFalse(V.shapes_are_node_local([broken])) + + @unittest.skipUnless(HAVE_PYSHACL, "pyshacl is in the image, not necessarily on the host") + def test_batched_and_whole_validation_agree(self): + """Same verdict and same violations, on a clean graph and a broken one.""" + from test import validation_fixtures as fixtures + + root = Path(vcf_rdfizer.__file__).parent + shapes, ontology = vcf_rdfizer.resolve_default_shacl_shapes(root, "core") + clean = fixtures.build_graph("expanded") + # Drop one record's POS: vcfc:pos is required, so this is a violation. + broken = "\n".join( + line for line in clean.splitlines() + if not (line.startswith("") and f"<{VCFC}pos>" in line) + ) + "\n" + for name, graph in (("clean", clean), ("broken", broken)): + with self.subTest(graph=name), tempfile.TemporaryDirectory() as td: + tmp = Path(td) + source = tmp / "g.nt" + source.write_text(graph, encoding="utf-8") + runs = {} + for label, rpb, workers in (("whole", None, 1), ("batched", 1, 2)): + out = tmp / label + out.mkdir() + runs[label] = V.validate_shacl( + source, shapes, out, ontology, records_per_batch=rpb, workers=workers) + whole, batched = runs["whole"], runs["batched"] + self.assertGreater(batched["batches"], 1) + self.assertEqual(batched["status"], whole["status"]) + self.assertEqual(batched["violationCount"], whole["violationCount"]) + self.assertEqual(sorted(batched["sample"]), sorted(whole["sample"])) + self.assertEqual(whole["status"], "PASS" if name == "clean" else "FAIL") + + +def _report_pid(_task): + """A shape task that only says which process ran it.""" + return True, f"pid {os.getpid()}\n" + + +def _report(*results): + """A pyshacl text report holding the given (severity, focus node) results.""" + blocks = "".join( + "Validation Result in MinCountConstraintComponent " + "(http://www.w3.org/ns/shacl#MinCountConstraintComponent):\n" + f"\tSeverity: sh:{severity}\n\tFocus Node: <{node}>\n" + "\tResult Path: vcfc:pos\n\tMessage: m\n" + for severity, node in results) + return f"Validation Report\nConforms: {not results}\nResults ({len(results)}):\n{blocks}" + + +class ShaclBatchMergeTests(VerboseTestCase): + """How batch outcomes become one report. The task is stubbed, so no pyshacl is needed.""" + + def setUp(self): + stub = mock.patch.dict(sys.modules, {"pyshacl": types.ModuleType("pyshacl")}) + stub.start() + self.addCleanup(stub.stop) + work = tempfile.TemporaryDirectory() + self.addCleanup(work.cleanup) + self.tmp = Path(work.name) + self.source = self.tmp / "g.nt" + self.source.write_text(ShaclBatchTests.GRAPH, encoding="utf-8") + + def validate(self, outcomes, source=None): + with mock.patch.object(V, "_validate_shacl_task", side_effect=outcomes) as task: + result = V.validate_shacl(source or self.source, [Path("shapes.ttl")], self.tmp, + records_per_batch=2, workers=1) + return result, task + + def test_a_result_repeated_across_batches_is_reported_once(self): + """Each batch carries the header, so its warnings recur in every batch.""" + header = _report(("Warning", "file://s.vcf#header/1")) + result, _task = self.validate([(False, header)] * 3) + self.assertEqual(result["batches"], 3) + self.assertEqual(result["advisoryCount"], 1) + self.assertEqual((result["status"], result["conforms"]), ("PASS", False), + "warnings alone do not fail the graph") + + def test_one_violating_batch_fails_the_graph(self): + bad = _report(("Violation", "file://s.vcf#record/3")) + result, _task = self.validate([(True, _report()), (False, bad), (True, _report())]) + self.assertEqual((result["status"], result["violationCount"]), ("FAIL", 1)) + self.assertIn("#record/3", result["sample"][0]) + + def test_a_batch_that_cannot_run_is_an_execution_failure(self): + result, _task = self.validate([(True, _report()), RuntimeError("parse failed"), (True, _report())]) + self.assertEqual(result["status"], "EXECUTION_FAILED") + self.assertIn("parse failed", result["error"]) + + def test_a_graph_without_records_is_validated_as_its_context(self): + header_only = self.tmp / "header.nt" + header_only.write_text(" .\n", encoding="utf-8") + result, task = self.validate([(True, _report())], source=header_only) + self.assertEqual(result["batches"], 1) + _context, batch, _shapes, _ontology = task.call_args.args[0] + self.assertIsNone(batch) + + def test_every_batch_gets_a_fresh_process(self): + """rdflib keeps freed memory, so a reused worker would hold its largest batch.""" + with mock.patch.object(V, "_validate_shacl_task", _report_pid): + result = V.validate_shacl(self.source, [Path("shapes.ttl")], self.tmp, + records_per_batch=1, workers=2) + pids = Path(result["report"]).read_text(encoding="utf-8").split()[1::2] + self.assertEqual((result["batches"], result["workers"]), (5, 2)) + self.assertEqual(len(set(pids)), 5) + self.assertNotIn(str(os.getpid()), pids) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/validation_fixtures.py b/test/validation_fixtures.py index 893cdf3..2c1ea2d 100644 --- a/test/validation_fixtures.py +++ b/test/validation_fixtures.py @@ -692,7 +692,7 @@ def parser_summary( digest[runner.record_digest_bucket([ record_iri, record.chrom, str(record.pos), record.record_id, record.ref, record.alt, - record.qual if include_qual else "", + runner.digest_qual(record.qual) if include_qual else "", record.filter_value, record.info, ])] += 1 summary["q11_record_digest"] = [ diff --git a/vcf_rdfizer.py b/vcf_rdfizer.py index 78bb7b8..0428fe1 100644 --- a/vcf_rdfizer.py +++ b/vcf_rdfizer.py @@ -4961,11 +4961,6 @@ def resolve_default_rules_path(repo_root: Path) -> Path: ) SHACL_PROFILE_CHOICES = ("core", "full") DEFAULT_SHACL_ONTOLOGY = "vcf-core-vocabulary.bundle.ttl" -#: pyshacl loads the whole graph into memory, so the default is size-gated -#: rather than unconditional. A fixture or a single-sample graph validates in -#: seconds; a cohort aggregate would not fit, and silently trying would turn a -#: safety net into an OOM. Above this, shapes stay available via --shacl-shapes. -DEFAULT_SHACL_MAX_SOURCE_BYTES = 512 * 1024 * 1024 def resolve_bundled_vocabulary_asset(repo_root: Path, relative: str) -> Path | None: @@ -5006,30 +5001,24 @@ def resolve_default_shacl_shapes( return shapes, ontology -#: The full profile set is quadratic-ish in graph size, so its gate is not the -#: core one. 16 MiB of VCF is a fixture or a small single-sample file, which is -#: where a two-minute structural check is a reasonable trade. +#: The full profile's SPARQL constraints compare records with each other, so it +#: is validated whole and grows quadratic-ish with the graph. 16 MiB of VCF is a +#: fixture or a small single-sample file, where a two-minute check is a +#: reasonable trade. The core profile has no such gate: the validator checks it +#: a batch of records at a time, so its memory does not follow the input. FULL_SHACL_MAX_SOURCE_BYTES = 16 * 1024 * 1024 -def shacl_max_source_bytes(profile: str) -> int: - """The size gate for one profile set.""" - return ( - FULL_SHACL_MAX_SOURCE_BYTES if profile == "full" - else DEFAULT_SHACL_MAX_SOURCE_BYTES - ) - - def shacl_default_applies(source_bytes: int | None, profile: str = "core") -> bool: """Whether to validate shapes by default for a source of this size. - Size-gated because pyshacl is in-memory. ``None`` means the size could not - be read, which is treated as too large: skipping a check is recoverable, - exhausting memory mid-run is not. + The core profile applies at any size. The full profile is size-gated, and + ``None`` -- an unreadable size -- counts as too large: skipping a check is + recoverable, exhausting memory mid-run is not. """ - if source_bytes is None: - return False - return 0 <= source_bytes <= shacl_max_source_bytes(profile) + if profile != "full": + return True + return source_bytes is not None and 0 <= source_bytes <= FULL_SHACL_MAX_SOURCE_BYTES def docker_image_exists(image: str) -> bool: @@ -5167,6 +5156,17 @@ def parse_positive_int(value: str, *, name: str) -> int: return parsed +def parse_non_negative_int(value: str, *, name: str) -> int: + """Parse a CLI integer where 0 is meaningful (it disables a limit).""" + try: + parsed = int(value) + except (TypeError, ValueError): + raise ValueError(f"{name} must be an integer") from None + if parsed < 0: + raise ValueError(f"{name} must be zero or a positive integer") + return parsed + + def compression_uses_partitioning(methods: list[str]) -> bool: """Return whether selected methods need bounded RDF chunks.""" return any(method in PARTITIONED_COMPRESSION_METHODS for method in methods) @@ -8976,6 +8976,8 @@ def run_validation_mode( # only place that knows the decoded graph's size. The size gate here # sees the packaged artifact and cannot. ("--shacl-max-triples", "shacl_max_triples"), + ("--shacl-batch-triples", "shacl_batch_triples"), + ("--shacl-workers", "shacl_workers"), ("--node-heap-mb", "node_heap_mb"), ): value = options.get(key) @@ -9527,8 +9529,8 @@ def main(): action="store_true", help=( "Skip the bundled SHACL shape layer, which is otherwise applied by " - "default to sources at or below " - f"{DEFAULT_SHACL_MAX_SOURCE_BYTES // (1024 * 1024)} MiB. The " + "default (core at any size, full up to " + f"{FULL_SHACL_MAX_SOURCE_BYTES // (1024 * 1024)} MiB). The " "default (core) profile checks structure the aggregate comparisons " "do not; --shacl-profile full adds the profiles that catch a " "corrupted value. --shacl-shapes overrides both" @@ -9589,11 +9591,25 @@ def main(): "--shacl-max-triples", default=None, help=( - "Skip the shape layer when the decoded graph exceeds this many " - "triples, recording the skip (0 disables). pyshacl is in-memory " - "and its cost tracks the graph, not the artifact it arrived in" + "Skip shapes validated whole when the decoded graph exceeds this " + "many triples, recording the skip (0 disables). pyshacl is " + "in-memory and its cost tracks the graph, not the artifact" ), ) + parser.add_argument( + "--shacl-batch-triples", + default=None, + help=( + "Validate node-level shapes in record batches of about this many " + "triples, so memory follows the batch, not the graph (0 validates " + "whole; default 500,000). SPARQL-based shapes are always whole" + ), + ) + parser.add_argument( + "--shacl-workers", + default=None, + help="Shape batches validated in parallel (default: up to 4)", + ) parser.add_argument( "--node-heap-mb", default=None, @@ -9765,16 +9781,19 @@ def main(): validation_engine_options["stop_after_query_timeout"] = True if args.validation_queries is not None: validation_engine_options["queries"] = args.validation_queries + # 0 is meaningful for both: no gate, and no batching. if args.shacl_max_triples is not None: - # 0 is meaningful here: it disables the gate, the way - # --validation-time-budget 0 means no ceiling. - try: - limit = int(args.shacl_max_triples) - except (TypeError, ValueError): - raise ValueError("--shacl-max-triples must be an integer") - if limit < 0: - raise ValueError("--shacl-max-triples must be zero or a positive integer") - validation_engine_options["shacl_max_triples"] = limit + validation_engine_options["shacl_max_triples"] = parse_non_negative_int( + args.shacl_max_triples, name="--shacl-max-triples" + ) + if args.shacl_batch_triples is not None: + validation_engine_options["shacl_batch_triples"] = parse_non_negative_int( + args.shacl_batch_triples, name="--shacl-batch-triples" + ) + if args.shacl_workers is not None: + validation_engine_options["shacl_workers"] = parse_positive_int( + args.shacl_workers, name="--shacl-workers" + ) if args.node_heap_mb is not None: validation_engine_options["node_heap_mb"] = parse_positive_int( args.node_heap_mb, name="--node-heap-mb" @@ -9822,8 +9841,8 @@ def main(): # a value that is counted but never read. Four of the ten mutation # classes the query suite misses are already covered by the # published profile, which no run in the benchmark campaign - # enabled. Default it on, size-gated, rather than leaving a - # written check permanently unused. + # enabled. Default it on (the full profile size-gated) rather than + # leaving a written check permanently unused. bundled_shapes, bundled_ontology = resolve_default_shacl_shapes( repo_root, args.shacl_profile )