diff --git a/.github/workflows/checks.yml b/.github/workflows/checks.yml new file mode 100644 index 0000000..7e43845 --- /dev/null +++ b/.github/workflows/checks.yml @@ -0,0 +1,81 @@ +# Checks that need no database. The query gates and the validators need a live +# TuringDB with the corpus loaded, which CI has no way to provide -- running +# them here would mean either skipping silently or building 85k structures on +# every push. They are listed in the README under Verification and are meant to +# be run against a live instance. +name: checks + +on: + push: + pull_request: + +jobs: + static: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: astral-sh/setup-uv@v3 + - run: uv sync --extra dev + + - name: ruff + run: uv run ruff check . + + - name: the CLI loads and its query library is well formed + run: | + uv run python - <<'PY' + import crystal + assert crystal.LIBRARY, "empty query library" + for name, (desc, cypher, tally) in crystal.LIBRARY.items(): + assert desc and cypher, name + assert cypher.upper().startswith(("MATCH", "CALL")), name + if tally: + for col in tally.split("+"): + assert col in cypher, f"{name}: tally column {col} not projected" + print(f"{len(crystal.LIBRARY)} named queries well formed") + PY + + - name: front-end destructure arity + run: | + uv run python - <<'PY' + import importlib.util, pathlib, sys + spec = importlib.util.spec_from_file_location( + "g", "tests/check_frontend_queries.py") + g = importlib.util.module_from_spec(spec); spec.loader.exec_module(g) + d = pathlib.Path("visualizer/src/components/viewer/crystal") + files = sorted(d.glob("*.ts")) + sorted(d.glob("*.tsx")) + bad = g.check_destructure_arity(files) + for b in bad: + print(b) + sys.exit(1 if bad else 0) + PY + + - name: the front-end query extractor still finds the shipped queries + run: | + uv run python - <<'PY' + import importlib.util, pathlib, sys + spec = importlib.util.spec_from_file_location( + "g", "tests/check_frontend_queries.py") + g = importlib.util.module_from_spec(spec); spec.loader.exec_module(g) + d = pathlib.Path("visualizer/src/components/viewer/crystal") + files = sorted(d.glob("*.ts")) + sorted(d.glob("*.tsx")) + q = g.extract(files) + missing = [t for t in g.REQUIRED_TOKENS + if not any(t in c for _, c in q)] + print(f"extracted {len(q)} queries") + if missing: + print("NOT COVERED: " + ", ".join(missing)) + sys.exit(1 if missing or len(q) < 40 else 0) + PY + + frontend: + runs-on: ubuntu-latest + defaults: + run: + working-directory: visualizer + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-node@v4 + with: + node-version: '20' + - run: npm ci + - run: npm run build diff --git a/.gitignore b/.gitignore index 8324960..e644c27 100644 --- a/.gitignore +++ b/.gitignore @@ -33,3 +33,9 @@ node_modules/ # ---- upstream visualizer bits that do not belong in this fork -------------- visualizer/.github/ visualizer/vercel.json + +# ---- internal working notes ------------------------------------------------ +# docs/DIALECT.md, docs/QUESTIONS.md and docs/PLAN_V2.md catalogue engine gaps +# and open chemistry judgement calls. Useful internally, wrong signal in a +# shared repo (see 059b8a9). Held on the build box instead. +docs/ diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..410ba90 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2026 TuringDB + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md index 9e287ca..11b5e05 100644 --- a/README.md +++ b/README.md @@ -1,145 +1,169 @@ # Crystallography on TuringDB -A demo built on the **Crystallography Open Database** (COD, CC0) that argues one -thing: crystal-structure data is natively a graph, and the most valuable layer -of it — the intermolecular contact network — is the layer nobody persists. +The intermolecular contact network of 84,805 crystal structures, stored as a +graph and queryable in milliseconds. COD and the CSD store atoms and coordinates. Hydrogen bonds, halogen bonds and molecular connectivity are *recomputed geometrically* every time somebody asks a question. This demo computes them once, at ingest, with the generating symmetry -operation attached to every edge, and then packing questions become pattern -matches instead of geometry jobs. +operation attached to every edge — so packing questions become pattern matches. -It also uses TuringDB's native commits: the corpus is ingested chronologically, -one commit per publication year, so `cod@2010` **is** the literature as it stood -at the end of 2010. A statistic published then still reproduces, without anyone -archiving a database dump. +Built from the Crystallography Open Database (CC0): CrystEngComm, Crystal Growth +& Design and Acta Cryst. E, the three journals closest to CCDC's small-molecule +remit. + +- **[TUTORIAL.md](TUTORIAL.md)** — the graph model, the queries, the UI, and how + the data gets in. Start there. --- ## Quick start ```bash -./run.sh # start TuringDB + the UI (data already built) -./run.sh --ingest # also acquire COD and build the graphs from scratch +./run.sh --ingest # first run: acquire COD and build the graphs +./run.sh # thereafter: start TuringDB + the UI on the built data ``` -Then open . +The CIFs are not in this repository — they are 2.8 GB and belong to COD — so the +first run has to be `--ingest`. It takes roughly 25 minutes. -Requirements: `uv`, Node 20+, and `rsync`. The `--ingest` path needs network -access to `crystallography.net`. Nothing else is assumed — the TuringDB server -binary ships inside the `turingdb` Python wheel, so there is no separate install. +Then open for the browser studio, or query from a shell: ---- +```bash +crystal facts # corpus summary +crystal run dimer # a named query +crystal query "MATCH (s:Structure) RETURN count(s)" # your own Cypher +crystal at 2009 "MATCH (s:Structure) RETURN count(s)" # as the corpus stood then +``` + +Requires `uv`, Node 20+ and `rsync`. The TuringDB server binary ships inside the +`turingdb` Python wheel, so there is nothing else to install. `--ingest` needs +network access to `crystallography.net`. -## What's in the graph +--- -Built from CrystEngComm, Crystal Growth & Design and Acta Cryst. E — the three -journals closest to CCDC's small-molecule remit. +## The corpus | | | |---|---| -| Structures | 9,899 | -| Atomic sites | 547,041 | -| Components (molecules & ions) | 10,766 | -| Covalent bonds | 584,812 | -| **Intermolecular contacts** | **33,930** | -| Space groups | 159 | +| Structures | 84,805 | +| Atomic sites | 4,694,045 | +| Components (molecules & ions) | 74,939 | +| Covalent bonds | 5,024,842 | +| **Intermolecular contacts** | **687,048** | +| Space groups | 214 | + +The contact layer is **three populations, never blended**, because a query has +to be able to refuse the weaker evidence: + +| | count | geometry | +|---|---|---| +| Strong H-bond, hydrogen located | 152,433 | mean H···A **2.048 Å** at **160.4°** | +| Inferred — no H refined, heavy-atom D···A fallback | 114,456 | 42.9% of the H-bond layer | +| Weak C–H···A, hydrogen located | 391,877 | mean H···A 2.674 Å at 145.0° | +| Halogen bond | 4,122 | C–X···A > 150°, X···A under ΣvdW | +| `close_contact` — excluded from every H-bond statistic | 24,160 | implausible inferred donor | + +82.1% of structures are in a centrosymmetric space group. + +Periodic-net dimensionality is scored for **40,320 structures (47.5%)** — the +rest have no located-hydrogen strong hydrogen bond and carry `net_dim = -1`. Any +dimensionality percentage quoted here is over that scored subset, not the whole +corpus. + +Every figure on this page is re-derived from the live graph by +`python tests/corpus_facts.py`. -Mean H···A distance **2.09 Å** at a mean D–H···A angle of **160.5°** — textbook -strong-hydrogen-bond geometry, which is the first sanity check a crystallographer -will run. +--- -### Graphs on the server +## What TuringDB is doing here that a relational store would not -| graph | what it is | -|---|---| -| `cod_slice_v2` | full schema: Structure, Atom, Component, Fragment, SpaceGroup, Element, Publication, Author, Journal | -| `cod_contacts_v2` | atoms + located-H hydrogen bonds only — the traversable projection (see *Known limits*) | -| `cod_versioned` | the same corpus committed chronologically, 26 commits, 2001→2026 | -| `cod_versions` | the tag ledger: `cod@` → commit hash | +**The contact network is persisted, not derived.** 687,048 edges, each carrying +the symmetry operation that generated the neighbour in CIF form (`2_565`). The +periodic network is therefore traversable without materialising a supercell, +and a motif search is a pattern match rather than a geometry job over every +candidate structure. ---- +**Molecular identity is cross-structure.** A `Component` node is deduplicated +corpus-wide by InChIKey, so a molecule appearing in forty structures is one node +with forty edges rather than forty rows joined on a string. Polymorph, solvate +and coformer questions are traversals from that node. Benzoic acid is one node; +so is water, shared by thousands of structures. + +**The motifs are cycles and paths.** An R²₂(8) acid dimer is a closed loop. +Periodic-net dimensionality is the rank of the lattice-translation subgroup +generated by cycles. Neither has a natural relational form; both are one query +here. + +**The corpus is versioned, not snapshotted.** The graph is ingested +chronologically, one commit per publication year. `cod@2009` *is* the +literature as it stood at the end of 2009 — the structures published later do +not exist in it. Reproducing a statistic published then is a checkout, not a +restore from an archived dump. -## The demo - -One studio page, reached from the flask icon in the left rail. - -**Left panel** — corpus overview, a colour legend for the nine node types, and -the full record of whatever node is selected on the canvas. - -**Investigate panel** (top right) — nine investigations. Each runs real Cypher, -paints the returned subgraph, reports nodes/edges/latency, exposes the Cypher, -and streams a narrative answer. Free text routes to the nearest investigation. - -Five of them carry a **commit toggle**: a strip of 26 ticks, one per publication -year, with bar height encoding corpus size. Click one and the same query re-runs -against the graph as it stood that year. - -### Suggested walkthrough - -1. **"How are space groups distributed in this corpus?"** — the opener. Purple - `SpaceGroup` hubs with blue `Structure` spokes, cluster labels naming the big - ones. Lands the textbook result: P2₁/c 34.9%, P-1 26.1%, centrosymmetric - groups 74.0%. -2. **Now click the commit ticks.** At `cod@2005` the corpus is 1,038 structures - and P2₁/c is 21.8%; at HEAD it is 9,899 and 34.9%. P2₁/c's dominance *grew* - over the period. This is the versioning story in one gesture — no archived - dumps, the corpus itself rewinds. -3. **"Which functional groups dominate the corpus?"** — the densest view, one - golden `Fragment` hub per group. 17,152 occurrences over 14 types. Mention - that 465 components carry a carboxylic acid and 565 a pyridine nitrogen: - that pair is the acid/pyridine heterosynthon of co-crystal design. -4. **Click a node** — the left panel reads the record back. A Structure gives - cell parameters, space group, Z, Z′, R factor and temperature; a Component - gives formula, InChIKey and solvent status. -5. **"Show the hydrogen-bond network of one structure"** — the contact layer - itself. Every edge carries the symmetry operation that generated the - neighbour, in CIF `2_565` form, which is what makes the periodic network - traversable without materialising a supercell. -6. **"If a structure is later corrected, which analyses are suspect?"** — the - compliance close. COD 2229029 enters at `cod@2011` and sits in 16 of 26 - commits: exactly the set of published analyses to re-examine. +Measured on this corpus (server execution time, median of three): + +| | | +|---|---| +| Count structures | 9 ms | +| Space-group filter | 59 ms | +| Packing subgraph, one structure | 65 ms | +| Three-hop synthon join | 118 ms | +| R²₂(8) dimer over 4.7 M atoms | 161 ms | +| `LOAD JSONL`, 21.2 M elements | 54.7 s (388 k elements/s) | +| Graph on disk | 2.85 GB — 32.8 kB per structure | --- -## Chemistry: what the criteria are +## Chemistry: the criteria Every cutoff lives in [`ingest/chemistry.py`](ingest/chemistry.py) with its literature source, and the SMARTS patterns in -[`ingest/fragments.py`](ingest/fragments.py), so a chemist can review both in -one place. - -- **Covalent bonds** — `d < r_cov(A) + r_cov(B) + 0.40 Å`, radii from - Cordero et al., *Dalton Trans.* 2008, 2832. Sites below 0.5 occupancy are - excluded so disordered partial sites are not bonded together. -- **Hydrogen bonds** — H···A < 2.5 Å with D–H···A > 120°, donors and acceptors - from {N, O, F, S, Cl}. Where no hydrogen was refined, a heavy-atom D···A - < 3.5 Å fallback is used and the edge is flagged `h_inferred`. The two - populations are never silently mixed. +[`ingest/fragments.py`](ingest/fragments.py), so both can be reviewed in one +place. + +- **Covalent bonds** — `d < r_cov(A) + r_cov(B) + 0.40 Å`, radii from Cordero + et al., *Dalton Trans.* 2008, 2832. Sites below 0.5 occupancy are excluded. +- **Hydrogen bonds** — H···A < 2.50 Å with D–H···A > 120°, donors and acceptors + from {N, O, F, S, Cl}. Chlorine appears as an acceptor almost entirely as + chloride: of 60 sampled Cl acceptor sites, none was covalently bonded to + carbon. +- **Inferred hydrogen bonds** — where no hydrogen was refined, a heavy-atom + D···A < 3.50 Å fallback applies and the edge is flagged `h_inferred`. This is + a **distance-only criterion: there is no angular test**, and 3.50 Å is loose + for O···O against the ~3.2 Å usually used. Only N, O and S are admitted as + inferred donors. The fallback fires only when the donor's *whole molecule* + carries no hydrogen — the "this entry has no H positions" case — rather than + per atom, which is what stops a carbonyl oxygen being treated as a donor in a + structure that is otherwise fully refined. 42.8% of the hydrogen-bond layer is + inferred; treat it as a separate, weaker population and filter it out of any + geometry claim. +- **Weak hydrogen bonds** — C–H···A with a *located* H, H···A ≤ 2.90 Å and + C–H···A ≥ 120°, after Taylor & Kennard, *JACS* **104** (1982) 5063 and + Steiner, *Angew. Chem. Int. Ed.* **41** (2002) 48. Stored as + `kind='hbond_weak'`; there is no inferred variant. - **Halogen bonds** — C–X···A with X ∈ {Cl, Br, I}, X···A below the sum of van der Waals radii (Bondi 1964, extended by Mantina 2009) and C–X···A > 150°. - **Components** — connected components of the covalent graph *after symmetry expansion*, so they are real molecules and ions rather than whatever the CIF - deposited. Identity is an InChIKey where RDKit perception succeeds, and a + deposited. Identity is an InChIKey where RDKit perception succeeds and a formula+connectivity digest otherwise, stored in a **separate property** so it is never mistaken for an InChIKey. ### Verification -Symmetry handling is the thing that would be spotted in seconds if it were -wrong, so it is tested rather than asserted: - ```bash -uv run python tests/validate_symmetry.py 200 # symop round-trip -uv run python tests/validate_chemistry.py 400 # valence plausibility -uv run python tests/validate_components.py 10 # eyeball component perception -uv run python tests/probe_turingdb_syntax.py # the Cypher surface +python tests/validate_symmetry.py 200 # symop round-trip +python tests/validate_chemistry.py 400 # valence plausibility +python tests/validate_components.py 10 # component perception +python tests/corpus_facts.py # re-derive every published figure +python tests/check_cli_queries.py # every CLI query runs +python tests/check_frontend_queries.py --negative-control # every UI query runs ``` `validate_symmetry` re-derives the operation stored on each edge and re-applies it: over 28,409 neighbour pairs, **100%** resolve and reproduce the neighbour's -position with a median residual of **1e-15 Å**. +position, median residual **1e-15 Å**. `validate_chemistry` reports valence distributions. In purely organic environments: H 98.2%, C 97.4%, N 99.6%, O 97.4%, S 97.5% plausible. @@ -151,97 +175,87 @@ has three bonds. ## Known limits -Stated up front, because they are the questions a crystallographer will ask. - - **π-stacking is not implemented.** It needs ring perception, and a - half-defined ring criterion would be worse than none. Hydrogen and halogen - bonds are in. -- **InChIKey perception reaches ~40% of components.** Metal complexes fail by - construction — InChI has no well-defined representation for them — and - components above 60 atoms are skipped because RDKit's bond-order search does - not terminate on them. Those keep the connectivity digest for identity. -- **Variable-length traversal cannot be restricted to an edge type.** TuringDB - rejects `-[e:CONTACT]->{1,8}` outright, so the contact network is materialised - as its own graph (`cod_contacts_v2`) where an untyped quantifier *is* a + centroid-distance criterion with no slippage constraint admits badly offset + pairs that most crystal engineers would not call stacking. +- **Hydrogen positions are not normalised.** X-ray H sits at ~0.85 Å (O–H) + against the neutron value of 0.983 Å, so the 2.50 Å H···A ceiling is + effectively ~2.35–2.40 Å normalised, and the measured H···A histogram is still + rising at that wall. D, H and A are all stored, so it can be renormalised + without re-running the geometry. +- **Disorder is an occupancy threshold, not disorder handling.** Sites below 0.5 + occupancy are excluded, but the commonest disorder ratio is exactly 50:50, so + both alternates survive. `_atom_site_disorder_group` is not read. +- **InChIKey perception reaches ~36% of components.** Metal complexes fail by + construction — InChI has no well-defined representation for them. +- **Halogen bonds are dominated by chlorine** (1,980 of 4,119 donors), and an + unactivated C–Cl is a very weak σ-hole donor. The 150° angle floor is loose + against the 155–165° convention — 458 edges, 11%, fall in that gap. Acceptors + are restricted to {N, O, F, S, Cl}, so X···π is undetectable; Cl···Cl **is** + detected and accounts for 1,184 edges. +- **`pyridine_nitrogen` misses the diazines.** The SMARTS requires the other five + ring atoms to be aromatic carbon, so pyrimidine, pyrazine, pyridazine and + triazine acceptors are excluded. `carboxylic_acid` requires an explicit + hydroxyl hydrogen, so acids in hydrogen-free refinements never match. +- **Z′ is COD metadata, not computed here.** +- **The `SUPERSEDES` layer is thin** — 684 edges, all asserted by COD's own + `duplicateof` column. Curated related-structure relationships would populate + it properly. +- **Variable-length traversal cannot be restricted to an edge type.** The engine + rejects `-[e:CONTACT]->{1,8}`, so the contact network is also materialised as + its own graph, `cod_contacts_v2`, where an untyped quantifier *is* a hydrogen-bond traversal. -- **Network dimensionality is a computation, not a traversal.** The stored - contact graph is a *quotient* graph, so a chain through the crystal is a short - cycle, not a long path — counting reachable nodes by depth reports everything - as isolated. `queries/dimensionality_net.py` does it properly, via the rank of - the lattice-translation subgroup generated by cycles. -- **The motif-drift numbers are not publishable from this corpus.** The - cumulative series is flat by construction and the per-year cohorts are too - thin at both tails (98 components in 2026). `queries/versioned_stats.py` - prints a do-not-quote warning below a 400-component threshold. The - *mechanism* is demonstrated; the finding would need the full 85k slice. -- **The narrative answers in the UI are not from a language model.** They are - templated prose over live `count()`/`avg()` queries. Every number is real and - re-derivable; nothing is invented. -- **Two space-group counts differ, deliberately.** `cod_slice_v2` keys a - `SpaceGroup` node by group *number* (159 nodes); `cod_versioned` keys by - Hermann-Mauguin *symbol* (195), so P2₁/c and P2₁/n are separate. Both are - correct and answer different questions. --- ## Provenance -CIFs come from the live rsync tree at `rsync://www.crystallography.net/cif/`, -not the published tarball — the newest release is `cod-rev297631-2025.02.09`, -roughly 19 months stale. Metadata comes from COD's public read-only MySQL -mirror. `data/manifest/slice_manifest.json` records the release revision, the -metadata SVN revision, the exact journal-name strings matched, the size cap -applied, and the sorted COD ID list with a checksum. +CIFs come from the live rsync tree at `rsync://www.crystallography.net/cif/`; +metadata from COD's public read-only MySQL mirror. COD's `robots.txt` is +`Disallow: /`, so nothing is scraped over HTTP. -Files above 512 kB are skipped (95.0% of the slice kept, 2.82 GB); the excluded -IDs are listed in `data/manifest/slice_skipped_ids.txt` rather than silently -dropped. COD's `robots.txt` is `Disallow: /`, so nothing is scraped over HTTP. +`--ingest` writes `data/manifest/slice_manifest.json` recording the release +revision, the metadata SVN revision, the exact journal-name strings matched, the +size cap applied, and the sorted COD ID list with a checksum; skipped IDs go to +`slice_skipped_ids.txt`. Those are **build outputs, not repository contents** — +`data/` is 2.8 GB and gitignored — so they describe your build, and two builds +from different COD revisions will differ. Files above 512 kB are skipped, which +keeps 95.0% of the slice. -COD data is dedicated to the public domain under CC0. +COD data is dedicated to the public domain under CC0. If you use COD in +published work, cite Gražulis et al., *J. Appl. Cryst.* **42** (2009) 726 and +*Nucleic Acids Res.* **40** (2012) D420. --- ## Layout ``` -ingest/ acquisition and the ingest pipeline +crystal.py the command-line client +ingest/ acquisition and the ingest pipeline download.py COD acquisition (rsync + MySQL) with an auditable manifest chemistry.py radii tables and geometric criteria, each cited fragments.py the SMARTS table - periodic.py explicit symmetry expansion + neighbour search - symmetry.py symop encoding in CIF form, and its verification + periodic.py symmetry expansion and neighbour search + symmetry.py symop encoding in CIF form structure.py bonds and periodic connected components - contacts.py hydrogen and halogen bond detection + contacts.py hydrogen, weak and halogen bond detection + netdim.py periodic-net dimensionality; the involution test perception.py RDKit identity and fragment matching + rdkit_guard.py killable worker for RDKit's bond-order search solvents.py solvent / counter-ion classification build_graph.py emits LOAD JSONL-shaped output - project_contacts.py the atoms+H-bonds projection + project_contacts.py the atoms+H-bonds projection load_turingdb.py load and verify - build_versioned.py one commit per publication year + the tag ledger -queries/ the analytical queries - dimensionality_net.py periodic-net dimensionality (0D/1D/2D/3D) - versioned_stats.py the same statistic across commits, plus the audit trail -tests/ validation and dialect probes -visualizer/ the front end (a trimmed TuringDB visualizer) - src/components/viewer/crystal/ the studio: panel, chat, canvas, labels + build_versioned.py one commit per publication year, plus the tag ledger +queries/ analytical scripts +tests/ validation and the two query gates +visualizer/ the browser studio ``` -Each ingest module's docstring records why it is shaped the way it is, including -the Cypher constructs this engine does and does not accept and the rewrite each -query needed. `tests/probe_turingdb_syntax.py` re-derives that surface against a -running instance if you want to check it yourself. - ---- - -## Credits - `visualizer/` is a trimmed fork of [turing-db/turingdb-visualizer](https://github.com/turing-db/turingdb-visualizer). -The WebGL canvas, node inspector, query bar and design system are theirs; this -repo adds `src/components/viewer/crystal/` and removes the studio pages that are -not about crystallography. +The WebGL canvas, node inspector and design system are theirs; this repo adds +`src/components/viewer/crystal/`. -Structural data is from the **Crystallography Open Database**, dedicated to the -public domain under CC0. If you use COD in published work, cite -Gražulis et al., *J. Appl. Cryst.* **42** (2009) 726 and -*Nucleic Acids Res.* **40** (2012) D420. +Licensed MIT — see [LICENSE](LICENSE). diff --git a/TUTORIAL.md b/TUTORIAL.md new file mode 100644 index 0000000..e47070f --- /dev/null +++ b/TUTORIAL.md @@ -0,0 +1,727 @@ +# Tutorial + +How the graph is built, what is in it, and how to ask it questions — from a +shell or from the browser. + +1. [Running it](#1-running-it) +2. [The graph model](#2-the-graph-model) +3. [Querying from the command line](#3-querying-from-the-command-line) +4. [The Cypher dialect: what works and what does not](#4-the-cypher-dialect-what-works-and-what-does-not) +5. [Worked questions](#5-worked-questions) +6. [The browser studio](#6-the-browser-studio) +7. [Time travel](#7-time-travel) +8. [How the data gets in](#8-how-the-data-gets-in) +9. [Writing fast queries](#9-writing-fast-queries) + +--- + +## 1. Running it + +```bash +git clone https://github.com/turing-db/crystallography +cd crystallography +./run.sh --ingest # first run only +``` + +The structures themselves are not in the repository: they are 2.8 GB of CIFs +belonging to COD, so the first run acquires them. It takes roughly 25 minutes +and needs network access to `crystallography.net`. Afterwards, `./run.sh` on its +own starts the server and the studio on the data already built. + +`run.sh` installs dependencies with `uv`, starts a TuringDB server on port 6691 +with the data directory `./turing-data`, builds the front end and serves it on +port 8087. The server binary ships inside the `turingdb` Python wheel; there is +no separate database to install. + +### Checking it is alive + +```bash +crystal graphs # graphs on the server, and whether each is loaded +crystal facts # corpus summary, every figure re-derived from the graph +``` + +If `crystal graphs` shows a graph as **not loaded**, load it before querying — +after a server restart, graphs are on disk but not resident, and an unloaded +graph answers `GRAPH_NOT_FOUND` rather than returning nothing: + +```python +from turingdb import TuringDB +c = TuringDB(host="http://localhost:6691") +c.load_graph("cod_slice_v2") # ~4 s for the full corpus +``` + +### The four graphs + +| graph | what it is | +|---|---| +| `cod_slice_v2` | the full model — this is the one you want | +| `cod_contacts_v2` | atoms + located-H hydrogen bonds only, for variable-length traversal | +| `cod_versioned` | the same corpus committed chronologically, one commit per publication year | +| `cod_versions` | the ledger mapping `cod@` to a commit hash | + +--- + +## 2. The graph model + +Nine node labels and twelve edge types, all populated. + +``` + ┌──────────────┐ + IN_SPACE_GROUP │ SpaceGroup │ hm_symbol, number, crystal_system, + ┌───────────────▶│ │ is_centrosymmetric + │ └──────────────┘ + │ ┌──────────────┐ + │ CONTAINS_ │ Element │ symbol, atomic_number + │ ELEMENT ────▶└──────────────┘ + │ +┌──┴────────────┐ HAS_SITE ┌──────────┐ IN_COMPONENT ┌─────────────┐ +│ Structure │────────────▶│ Atom │───────────────▶│ Component │ +│ │ │ │ │ │ +│ cod_id │ │ uid │ │ inchikey │ +│ hm_symbol │ │ label │ │ formula │ +│ a b c α β γ │ │ element │ │ charge │ +│ Z, Z_prime │ │ fract_* │ │ is_solvent │ +│ r_factor │ │ occupancy│ │ is_polymeric│ +│ temperature │ └──────────┘ └─────────────┘ +│ net_dim │ │ │ │ +│ net_dim_weak │ BONDED_TO CONTACT HAS_FRAGMENT +│ has_hydrogens │ │ │ ▼ +└───────────────┘ ▼ ▼ ┌─────────────┐ + │ │ (Atom) (Atom) │ Fragment │ + │ │ │ fragment_ │ + │ └── CONTAINS_COMPONENT ──▶ (Component) │ type │ + │ role: principal | solvent | │ smarts │ + │ counter_ion └─────────────┘ + │ + ├── PUBLISHED_IN ──▶ (Publication) ── IN_JOURNAL ──▶ (Journal) + │ │ + │ └── AUTHORED_BY ──▶ (Author) + └── SUPERSEDES ────▶ (Structure) +``` + +### The edge that matters: `CONTACT` + +Every intermolecular contact is an edge between two `Atom` nodes, carrying: + +| property | meaning | +|---|---| +| `kind` | `hbond` · `hbond_weak` · `halogen` · `close_contact` | +| `h_inferred` | `true` when no hydrogen was refined and a heavy-atom cutoff was used | +| `length` | H···A, or X···A for a halogen bond, or D···A when `h_inferred` | +| `angle` | D–H···A, or C–X···A | +| `symop` | the generating operation in CIF form, e.g. `2_565` | +| `symop_triplet` | the same operation written out, e.g. `-x+1,-y+1,-z+1` | +| `is_involution` | `true` if the operation is its own inverse | +| `cod_id` | the structure the contact belongs to | + +**Always filter on `kind` and `h_inferred`.** The three populations are +deliberately separate and blending them silently is how packing statistics go +wrong. Nearly every query in this tutorial begins: + +```cypher +WHERE r.kind = 'hbond' AND r.h_inferred = false +``` + +**`is_involution` separates a ring motif from a chain.** A centrosymmetric dimer +is generated by an inversion, mirror or two-fold — an operation that is its own +inverse, so it relates exactly two molecules. A catemer is generated by a 2₁ +screw or a glide, which has infinite order and builds a chain. Both have a +two-fold *rotation* part, so nothing short of composing the operation with +itself distinguishes them. It is precomputed at ingest. + +**The identity operation passes this test**, and about 17% of the hits in the +dimer query carry `symop = '1_555'`. Those are contacts between two +crystallographically independent molecules in a Z′ > 1 cell — real pairs, but +generated by no symmetry at all. Add `AND h.symop <> '1_555'` if you want only +the symmetry-generated motif, and be aware that doing so discards genuine +general-position dimers. + +### Two things that are easy to get wrong + +**`Fragment` is a 16-node lookup table, and it is reached two ways.** +`HAS_FRAGMENT` runs Component → Fragment and means *this molecule contains at +least one of these somewhere*. `IN_FRAGMENT` runs Atom → Fragment and means +*this specific atom is part of one*. + +Use the right one. A motif query must use `IN_FRAGMENT`: written through +`HAS_FRAGMENT`, `f.fragment_type = 'carboxylic_acid' AND a1.element = 'O'` +means "some oxygen of a molecule that has a COOH somewhere", which in a +polyfunctional molecule also matches its alcohols, esters and nitro oxygens. +A question genuinely about molecules — "what co-crystallises with an acid" — +should use `HAS_FRAGMENT`. + +**`Component` is deduplicated across structures.** One node per distinct +molecular identity for the whole corpus, so `c.formula` is that shared node's +formula and `CONTAINS_COMPONENT` fans out to every structure containing it. +That is what makes cross-structure questions one hop, and it is also why any +pattern reaching two atoms through two Component nodes needs scoping to a +structure. + +### Two things about the representation + +**The stored graph is a quotient graph.** Nodes are asymmetric-unit sites, and +each contact edge carries the operation to its neighbour. A hydrogen-bonded +chain running through the crystal is therefore a short *cycle*, not a long path. +Counting reachable nodes by depth reports everything as isolated; see +`ingest/netdim.py` for the correct treatment. + +**A centrosymmetric dimer is one edge, not two.** Its second hydrogen bond is +the symmetry image of the first, so the quotient graph stores it once. If you +are counting physical hydrogen bonds rather than stored edges, an +involution-generated edge counts twice. + +--- + +## 3. Querying from the command line + +```bash +crystal facts # corpus summary +crystal list # the named queries +crystal run dimer # run one +crystal run dimer --cypher # print its Cypher instead of running it +crystal query "MATCH (s:Structure) RETURN count(s)" +crystal at 2009 "MATCH (s:Structure) RETURN count(s)" +crystal commits # the version ledger +crystal graphs # graphs on the server +``` + +Output options work before or after the subcommand: + +```bash +crystal run dimer --limit 0 --csv > dimers.csv +crystal query "MATCH (s:Structure) RETURN s.cod_id, s.r_factor" --json +crystal --graph cod_contacts_v2 query "MATCH (n:Atom) RETURN count(n)" +``` + +`TURING_HOST` and `TURING_GRAPH` set the defaults. + +The named queries are a starting point, not a fixed menu — `--cypher` prints the +query so you can paste it into `crystal query`, change a predicate and re-run. + +### From Python + +```python +from turingdb import TuringDB + +c = TuringDB(host="http://localhost:6691") +c.load_graph("cod_slice_v2") # only needed once per server start +c.set_graph("cod_slice_v2") + +df = c.query("MATCH (s:Structure) WHERE s.net_dim = 2 RETURN count(s)") +print(df) # results come back as a DataFrame +``` + +### Over REST + +Useful from any language. The body is the raw Cypher string, not JSON-wrapped. + +```bash +curl -s -X POST "http://localhost:6691/query?graph=cod_slice_v2" \ + -H "Content-Type: application/json" \ + -d "MATCH (s:Structure) RETURN count(s)" +``` + +The response is **column-oriented and chunked**: + +```json +{"header": {"column_names": ["count(s)"], "column_types": ["UInt64"]}, + "data": [[[84801]]], "time": 0.144} +``` + +`data` is a list of chunks; each chunk is a list of *columns*, not rows. To +rebuild rows from chunk `k`: read `n = chunk[0].length`, then +`row[i][c] = chunk[c][i]`. + +**Failures arrive as HTTP 200 with a non-null top-level `error`.** Checking the +status code is not enough. + +--- + +## 4. The Cypher dialect: what works and what does not + +TuringDB implements most of Cypher. The gaps that matter in practice: + +**Not supported.** `WITH`, `collect()`, `OPTIONAL MATCH`, `UNION`, `DISTINCT`, +`MERGE`, variable-length path aliases, `IN`, `CONTAINS` / `STARTS WITH` / `=~`. + +**Only two aggregates exist.** `count()` and `avg()` work; `min()`, `max()`, +`sum()` and `stdev()` are all rejected. So "the shortest contact", "the longest +cell axis" and "the total number of sites" cannot be expressed as aggregates — +project the column and reduce it client-side, which is what `crystal` does. + +**No GROUP BY, and aggregates do not combine with other return items.** +`RETURN f.fragment_type, count(c)` is rejected, and so is an aggregate inside +`ORDER BY`. A "how many of each" question is expressed as a flat projection and +counted client-side — which is exactly what `crystal run fragments` does, and +why `--cypher` shows a query with no `count()` in it. + +**`IS NULL` does not work.** Every optional numeric field has an explicit +boolean companion instead: `has_temperature`, `has_r_factor`, `has_net_dim`, +`has_inchikey`. Absent properties still project safely as nulls, and comparisons +silently exclude them. + +**There are no query parameters.** `$name` is a `PARSE_ERROR: Not implemented: +Parameters`, so a value can only reach a query by string interpolation — which +means **the usual protection against injection is not available**: + +```cypher +WHERE g.hm_symbol = 'P -1' -- 1 +WHERE g.hm_symbol = 'P -1' OR 1=1 -- 214 +``` + +Every query in this repository interpolates, because there is no alternative. +That is acceptable here: the corpus is public CC0 data, the server is bound to +localhost, and the demo has no untrusted input. **It is not acceptable in +anything that accepts input from a user**, and the patterns in this repository +should not be lifted into such a service without an escaping or allow-listing +layer in front of them. Until parameters land, treat a Cypher string as +something you build, never something a caller supplies. + +**`type` is a reserved word.** The edge-type function is `edgeType(r)`, and a +property named `type` needs backticks. This is why fragments carry +`fragment_type`. + +**Property types are global, not per label.** One property name has one type +across the whole graph. `Publication.volume` is `journal_volume` here because +`Structure.cell_volume` already claimed a numeric `volume`. + +**Projecting a property that no node carries is an `ANALYZE_ERROR`**, not an +empty column. Same for an edge type that does not exist at the commit you are +querying. Use `CALL db.propertyTypes()` and `CALL db.edgeTypes()` to see what +exists — but note they read a registry that remembers names from earlier builds, +so they over-report; confirm with a `count()` before relying on one. + +**Variable-length traversal cannot be restricted to an edge type.** The +quantifier is postfix — `-[e]->{1,8}`, `-[e]->+`, `-[e]->*` — and adding a type +is rejected. That is why `cod_contacts_v2` exists: in a graph whose only edges +are hydrogen bonds, an untyped quantifier *is* a hydrogen-bond traversal. + +```cypher +MATCH (a:Atom {uid:'2229029_O1'})-[e]->{1,6}(b:Atom) RETURN count(b) +``` + +**`IN` does not work, but `UNWIND` does** — and it is the replacement: + +```cypher +UNWIND [1, 2] AS k MATCH (s:Structure) WHERE s.net_dim = k RETURN count(s) +``` + +**Useful and easy to miss:** `SKIP`, numeric comparisons in `WHERE`, arithmetic +in `RETURN` (`RETURN r.length * 2`), `labels(n)`, `count(*)`, multi-key +`ORDER BY`, `ORDER BY `, and `cosine_similarity()` / +`euclidean_distance()` over embedding properties. + +**`toInteger()` and `toFloat()` do not work** despite appearing in the Cypher +surface — both are an `ANALYZE_ERROR` here. + +--- + +## 5. Worked questions + +### Which space groups dominate? + +No GROUP BY, so project flat and tally: + +```bash +crystal run spacegroups +``` + +```cypher +MATCH (s:Structure)-[:IN_SPACE_GROUP]->(g:SpaceGroup) RETURN g.hm_symbol +``` + +`hm_symbol` exists on **both** labels and means different things, which is the +single easiest way to get a wrong answer here. + +`SpaceGroup` is keyed by group **number** — one node per group — and its +`hm_symbol` is whichever setting happened to be seen first. So this returns 0, +because the node for No. 14 carries `P 1 21/n 1`: + +```cypher +MATCH (s:Structure)-[:IN_SPACE_GROUP]->(g:SpaceGroup) +WHERE g.hm_symbol = 'P 1 21/c 1' RETURN count(s) -- 0, and not what you meant +``` + +Ask by **number** when you want the group, and use `Structure.hm_symbol` when +you want the setting as the depositor reported it: + +```cypher +MATCH (s:Structure)-[:IN_SPACE_GROUP]->(g:SpaceGroup) +WHERE g.number = 14 RETURN count(s) -- 29,644, all settings + +MATCH (s:Structure) WHERE s.hm_symbol = 'P 1 21/c 1' +RETURN count(s) -- 18,643, that setting +``` + +For the same reason, `crystal run spacegroups` tallies by group and labels each +bucket with one representative symbol: the 29,644 printed against `P 1 21/n 1` +is all of No. 14, not that setting alone. `cod_versioned` keys `SpaceGroup` by +symbol instead, so there the settings are separate nodes. + +### The hydrogen-bond network of one structure + +```cypher +MATCH (a:Atom)-[r:CONTACT]->(b:Atom) +WHERE a.cod_id = 2229029 AND r.kind = 'hbond' +RETURN a.label, b.label, r.length, r.angle, r.symop_triplet +``` + +Every row carries the operation that generated the neighbour, so the periodic +network is reconstructible without building a supercell. + +### R²₂(8) carboxylic-acid dimers + +The canonical graph-set motif. Walk molecule → atom → contact → atom → molecule, +and require the generating operation to be an involution: + +```cypher +MATCH (f1:Fragment)<-[:HAS_FRAGMENT]-(c1:Component)<-[:IN_COMPONENT]-(a1:Atom) + -[h:CONTACT]->(a2:Atom)-[:IN_COMPONENT]->(c2:Component) + -[:HAS_FRAGMENT]->(f2:Fragment) +WHERE f1.fragment_type = 'carboxylic_acid' + AND f2.fragment_type = 'carboxylic_acid' + AND h.kind = 'hbond' AND h.h_inferred = false + AND h.is_involution = true + AND a1.element = 'O' AND a2.element = 'O' + AND a1.cod_id = a2.cod_id +RETURN a1.cod_id, a1.label, a2.label, h.symop_triplet, h.length, h.angle +``` + +Swap `is_involution = true` for `false` and the same query returns C(4) +catemers instead — screw axes and glides rather than inversion centres. + +`a1.cod_id = a2.cod_id` is redundant — `a1` and `a2` are bound by the same +`CONTACT` edge, which never crosses a structure — but it is kept as an explicit +reminder that `Component` nodes *are* deduplicated corpus-wide, so any pattern +that reaches two atoms through two Component nodes rather than through one edge +does need scoping. + +### Acid versus pyridine: which synthon wins? + +```bash +crystal run synthon-homo # O-H...O acid to acid +crystal run synthon-hetero # O-H...N acid to pyridine +``` + +Both counts are stored *edges*. To compare physical hydrogen bonds, count +involution-generated homosynthon edges twice — a centrosymmetric acid dimer is +two bonds held as one edge, and the heterosynthon has no such pairing. On this +corpus that moves the split from **62.2 / 37.8** to **75.0 / 25.0**: the +homodimer wins either way, and the correction widens the margin rather than +changing the answer. + +Note both figures are over a corpus where only 110 molecules carry *both* +groups, so this is the corpus-wide propensity of each synthon, not a +head-to-head contest within molecules that had a genuine choice. Scope it to +those 110 if that is the question you mean. + +### Is a packing a dimer, a chain, a sheet or a framework? + +Precomputed at ingest as `net_dim` (strong bonds) and `net_dim_weak` (weak +admitted as well): `0` finite motif, `1` chain, `2` sheet, `3` framework. + +`-1` means **not scored** — the structure has no located-hydrogen strong +hydrogen bond, or the orbit expansion exceeded its ceiling. That is **44,486 of +84,805 structures, 52.5% of the corpus**, so any dimensionality percentage is +over the 40,320 that were scored, not over the whole corpus. Filter on +`has_net_dim = true` to be explicit about it. (`-2` is reserved for a +computation that raised; it is currently empty.) + +```cypher +MATCH (s:Structure) WHERE s.net_dim = 2 RETURN count(s) +``` + +`net_dim_is_framework` marks structures whose *covalent* net is already periodic +— coordination polymers, where the number describes the framework rather than +the hydrogen-bond net. Exclude them if that matters to your question. + +### Where else does this molecule appear? + +```cypher +MATCH (c:Component)<-[:CONTAINS_COMPONENT]-(s:Structure) +WHERE c.inchikey = 'WPYMKLBDIGXBTP-UHFFFAOYSA-N' +RETURN s.cod_id, s.hm_symbol, s.year +``` + +That one returns **82 structures across 13 different space-group settings** — +one node, 82 edges, and the polymorphism falls out of the same traversal. In a +relational schema this is a self-join over a structure table on a chemical +identity string, which only works if that string was computed consistently in +the first place, and still says nothing about the packing without a second pass. + +`crystal run recurrence` ranks molecules by how many structures they appear in. + +--- + +## 6. The browser studio + +, then pick a crystallography graph. + +**Left panel** — corpus overview with the three contact populations separated, a +colour legend for the nine node types, and the full record of whatever node is +selected on the canvas. + +**Investigate panel** (top right) — a library of queries on three shelves: + +- *Packing & motifs* — the contact network itself +- *Solid form & design* — salt vs co-crystal, coformers, polymorphs +- *Provenance & versioning* — the commit history + +Each entry says what it demonstrates, runs real Cypher, paints the returned +subgraph, and reports nodes, edges and latency. **show Cypher** exposes the +query; **edit & run this query** drops it into the input box so you can change a +predicate and see the canvas move. + +**The input box takes Cypher directly.** Anything beginning `MATCH`, `CALL`, +`RETURN` or similar is executed; everything else routes to the nearest canned +investigation. A query projecting the eight-column subgraph shape + +``` +RETURN a, labels(a), , e, edgeType(e), b, labels(b), +``` + +paints the canvas. Anything flatter — a count, an average, a projection — comes +back as a table. Engine errors are shown verbatim, because `PARSE_ERROR`, +`ANALYZE_ERROR` and `PLAN_ERROR` mean quite different things. + +The standard visualizer toolbar is hidden on these graphs on purpose: its data +pipeline fetches every node and edge in the graph, which is not viable at 4.7 M +atoms, and it would overwrite the subgraph the studio paints. + +--- + +## 7. Time travel + +The corpus is ingested chronologically, one commit per publication year: 26 +years, and 26 entries in the `cod_versions` ledger. + +```bash +crystal commits # year, tag, totals, commit hash +crystal at 2009 "MATCH (s:Structure) RETURN count(s)" # 32,419 +crystal at 2026 "MATCH (s:Structure) RETURN count(s)" # 84,805 +``` + +In the studio, five investigations carry a strip of 26 ticks below the answer; +clicking one re-runs the same query against that commit. + +Over REST, two steps — a commit must be made resident before it can be read: + +```bash +curl -s -X POST "http://localhost:6691/query?graph=cod_versioned" \ + -H "Content-Type: text/plain" -d "LOAD COMMIT ''" +curl -s -X POST "http://localhost:6691/query?graph=cod_versioned&commit=" \ + -H "Content-Type: text/plain" -d "MATCH (s:Structure) RETURN count(s)" +``` + +`LOAD COMMIT` is idempotent. + +Two properties of the versioning worth knowing: + +- **A commit survives a restart; an open change does not.** Anything built on an + unsubmitted change has to be re-created, and there is no way to list changes. +- **Commit hashes are build-time values.** Rebuilding the same graph from the + same script produces different hashes, so never hardcode one. Resolve it at + runtime from the `cod_versions` ledger. +- **`CALL db.history()` returns 54 entries, not 26.** Each publication year + costs two commits — an intermediate `COMMIT` so that the year's edge writes + can see the nodes they attach to, then the `CHANGE SUBMIT` — plus the initial + one. The 26 in the ledger are the ones worth checking out; use the ledger + rather than picking entries out of the history. + +**What is versioned.** `cod_versioned` carries `Structure`, `SpaceGroup`, +`Component` and `Fragment`, joined by `IN_SPACE_GROUP`, `CONTAINS_COMPONENT` and +`HAS_FRAGMENT`. It deliberately holds **no atoms and no contacts** — that layer +is 5.4 M edges and committing it 26 times would multiply the corpus by an order +of magnitude for no gain in the question it answers. + +So structure counts, space-group distributions and fragment censuses rewind; +hydrogen-bond geometry, dimer counts and dimensionality do not. Asking for them +at a commit is an honest error rather than a wrong number: + +``` +$ crystal at 2009 "MATCH ()-[r:CONTACT]->() RETURN count(r)" +error: ANALYZE_ERROR: Unknown edge type: CONTACT +``` + +TuringDB has commit hashes, integer change ids and a single `main` — no tags and +no branches. The `cod@` names live in the `cod_versions` ledger graph, +which makes the mapping itself queryable. + +--- + +## 8. How the data gets in + +``` +COD rsync + MySQL → parse cache → JSONL → LOAD JSONL → cod_slice_v2 + download.py build_graph.py load_turingdb.py +``` + +### Acquisition + +`ingest/download.py` fetches CIFs from `rsync://www.crystallography.net/cif/` +with `--files-from`, so it pulls the ~2 GB slice rather than the 109 GB tree, +and metadata from COD's public read-only MySQL mirror. It writes a manifest +recording the release revision, the SVN revision at download time, the exact +journal strings matched, the size cap, and the sorted ID list with a checksum. + +### Parsing + +`ingest/build_graph.py` parses each CIF with gemmi, then per structure: + +1. builds atoms from the asymmetric unit; +2. expands symmetry and searches for neighbours (`periodic.py`); +3. detects covalent bonds and finds connected components — real molecules and + ions, after expansion (`structure.py`); +4. detects hydrogen, weak and halogen contacts (`contacts.py`); +5. perceives identity and functional groups with RDKit (`perception.py`); +6. computes periodic-net dimensionality (`netdim.py`). + +Results are cached per structure, so `--from-cache` re-emits the graph without +re-parsing. Bond-order perception runs in a killable child process +(`rdkit_guard.py`) because RDKit's search does not always terminate and is a C++ +call no Python watchdog can interrupt. + +### Loading + +The graph is written as JSONL in Neo4j APOC's `apoc.export.json.all` shape — one +object per line, nodes and relationships interleaved: + +```json +{"type":"node","id":"0","labels":["Structure"],"properties":{"cod_id":1100119}} +{"type":"relationship","id":"0","label":"HAS_SITE","start":{"id":"0"},"end":{"id":"1"},"properties":{}} +``` + +**Node ids and relationship ids must each start at 0 and increment with no +gaps.** A crc32-style identifier has to be renumbered into a dense index at +export; keep the original as an ordinary property. + +The file goes in `/data/`, then: + +```cypher +LOAD JSONL 'cod_slice_v2.jsonl' AS cod_slice_v2 +``` + +Two practical notes. `LOAD JSONL` **only ever creates** — loading onto an +existing graph is an `EXEC_ERROR`, and there is no `DROP GRAPH`, so rebuilding +means stopping the server and removing `/graphs/`. And the +statement must be issued with `?graph=` pointing at an *existing* graph; +targeting the name being created is `GRAPH_NOT_FOUND`. + +Serialise with `allow_nan=False`. Python's `json.dumps` emits a bare `NaN` by +default, which is not valid JSON, and one occurrence rejects the entire file — +the API reports only `Failed to load JSONL graph`, with the line and column in +`/logs/`. + +### The two write paths, and how to choose + +| | `LOAD JSONL` | Cypher `CHANGE NEW … COMMIT … CHANGE SUBMIT` | +|---|---|---| +| Same 84,805 structures | **55 s** | **77 min** | +| History | none — one immutable graph | 26 queryable commits | +| Cost profile | linear | quadratic in corpus size | + +The chronological build is slower because each year's writes `MATCH` against a +graph that grows, and the dialect has no property index on that graph. Pick by +whether you need the history, not by preference. `ingest/build_versioned.py` +does the second. + +### Property indexes + +```cypher +CREATE INDEX ix_frag_type FOR (n) ON n.fragment_type +CREATE INDEX ix_contact_kind FOR [e] ON e.kind +``` + +These are **writes** — they must run inside a change. `CALL db.showIndexes()` +lists them. The optimiser rewrites `WHERE n.prop = ` into an index lookup +automatically. They are worth having, but see the next section: they are a +tiebreaker, not the main lever. + +--- + +## 9. Writing fast queries + +**Write a multi-hop read as one connected chain, starting at the most selective +label.** This is worth more than everything else combined. + +Slow — separate patterns joined in the `WHERE`, which forces a cross product: + +```cypher +MATCH (a1:Atom)-[h:CONTACT]->(a2:Atom), + (a1)-[:IN_COMPONENT]->(c1:Component)-[:HAS_FRAGMENT]->(f1:Fragment), + (a2)-[:IN_COMPONENT]->(c2:Component)-[:HAS_FRAGMENT]->(f2:Fragment) +WHERE f1.fragment_type = 'carboxylic_acid' AND f2.fragment_type = 'carboxylic_acid' +``` + +Fast — one chain, beginning at `Fragment`, of which there are 16: + +```cypher +MATCH (f1:Fragment)<-[:HAS_FRAGMENT]-(c1:Component)<-[:IN_COMPONENT]-(a1:Atom) + -[h:CONTACT]->(a2:Atom)-[:IN_COMPONENT]->(c2:Component) + -[:HAS_FRAGMENT]->(f2:Fragment) +WHERE f1.fragment_type = 'carboxylic_acid' AND f2.fragment_type = 'carboxylic_acid' +``` + +Same rows. On this corpus the first does not return at all — it was stopped at +900 s — and the second runs in **157 ms**. At one-eighth the corpus size the +same pair is 17.2 s against 0.84 s, so the penalty grows faster than the data. +Property indexes on the slow form bought 30%; the rewrite is worth four orders +of magnitude. + +The planner follows the pattern as written. Start from whichever label has the +fewest nodes and chain outward. + +**Other things that help:** + +- Filter on `kind` early: it is indexed (`CALL db.showIndexes()` lists what is). + `h_inferred` is not indexed, but it is cheap once `kind` has cut the set. +- Prefer the precomputed property to recomputing: `net_dim` and `is_involution` + exist so that dimensionality and ring-vs-chain cost nothing at query time. +- Aggregates over large edge sets are the slowest common operation; `avg()` over + 152,592 edges takes ~1.1 s where every traversal here is under 250 ms. +- For exploration, use `CALL db.labels()`, `CALL db.propertyTypes()` and + `CALL db.edgeTypes()` rather than `MATCH (n) RETURN n`. + +### `shortestPath` + +It works, and for a contact network it is the obvious thing to reach for: + +```cypher +MATCH (a:Atom), (b:Atom) shortestPath(a, b, w, dist, path) RETURN dist +``` + +**Do not project an endpoint variable.** `RETURN a.label` — a String property of +an endpoint — terminates the server process. `RETURN a` and `RETURN a.cod_id` +are a clean `PLAN_ERROR`; only `RETURN dist` and `RETURN dist, path` are safe. +`path` is a positionally interleaved list of node and edge ids whose two id +spaces overlap, so resolve the even positions with `CALL db.getNodes`. + +Given that `label` is a String on every `Atom` here, this is easy to hit. Treat +it as a sharp edge rather than a feature until it is fixed upstream. + +--- + +## Checking your changes + +Two gates run the queries the project ships against a live server, and both +fail loudly rather than returning empty: + +```bash +python tests/check_cli_queries.py +python tests/check_frontend_queries.py --negative-control +``` + +The front-end gate recovers queries by scanning the source, so it also asserts +that a set of required schema elements — `CONTACT`, `IN_COMPONENT`, +`HAS_FRAGMENT`, `net_dim`, `is_involution`, `Snapshot`, `IN_SPACE_GROUP` — +appears somewhere in what it extracted, and fails if any is missing. Without +that, a scraper that quietly stops matching reports success over an empty set. + +It also checks that every `const [...] = await Promise.all([...])` names exactly +as many variables as the array has entries: the panels read their results +positionally, so a mismatch renders a plausible wrong number rather than an +error. + +```bash +python tests/corpus_facts.py # re-derive every figure the README quotes +``` diff --git a/crystal.py b/crystal.py new file mode 100755 index 0000000..a45fecc --- /dev/null +++ b/crystal.py @@ -0,0 +1,515 @@ +#!/usr/bin/env python3 +"""Command-line access to the crystallography graph. + +Everything the browser studio can ask, this can ask too. The studio is better +for looking at a packing motif; this is better for a corpus-wide number you +want to paste into a paper, a loop over many structures, or a quick check that +a criterion does what you think it does. + + crystal facts corpus summary, all figures re-derived + crystal list the named queries + crystal run dimer run one by name + crystal run dimer --cypher print its Cypher without running it + crystal query "MATCH (s:Structure) RETURN count(s)" + crystal query "..." --csv > out.csv + crystal at 2009 "MATCH (s:Structure) RETURN count(s)" + crystal commits the version ledger + +Reads TURING_HOST (default http://localhost:6691) and TURING_GRAPH +(default cod_slice_v2). +""" + +from __future__ import annotations + +import argparse +import csv +import json +import os +import sys +import urllib.error +import urllib.parse +import urllib.request + +DEFAULT_HOST = os.environ.get("TURING_HOST", "http://localhost:6691") +DEFAULT_GRAPH = os.environ.get("TURING_GRAPH", "cod_slice_v2") +VERSIONED_GRAPH = "cod_versioned" +LEDGER_GRAPH = "cod_versions" + + +# -------------------------------------------------------------------------- +# transport +# -------------------------------------------------------------------------- + + +class QueryError(RuntimeError): + pass + + +def query(cypher: str, graph: str = DEFAULT_GRAPH, host: str = DEFAULT_HOST, + commit: str | None = None) -> tuple[list[str], list[list]]: + """Run one Cypher statement. Returns (column names, rows). + + TuringDB answers column-oriented and in chunks: the payload is + `{header: {column_names, column_types}, data: [chunk, ...]}` where each + chunk is a list of COLUMNS, not rows. Failures arrive with HTTP 200 and a + non-null top-level `error`, so the status code alone is not a check. + """ + url = f"{host}/query?graph={urllib.parse.quote(graph)}" + if commit: + url += f"&commit={urllib.parse.quote(commit)}" + req = urllib.request.Request( + url, data=cypher.encode(), + headers={"Content-Type": "application/json"}, method="POST") + try: + with urllib.request.urlopen(req, timeout=900) as r: + payload = json.loads(r.read()) + except urllib.error.URLError as exc: + raise QueryError( + f"cannot reach TuringDB at {host} ({exc.reason}). " + f"Is the server running? See README, 'Running it'.") from None + if payload.get("error"): + detail = payload.get("error_details") or "" + raise QueryError(f"{payload['error']}: {detail}".strip().rstrip(":")) + + names: list[str] = payload.get("header", {}).get("column_names", []) + cols: list[list] = [] + for chunk in payload.get("data") or []: + if not isinstance(chunk, list): + continue + for i, col in enumerate(chunk): + if not isinstance(col, list): + continue + while len(cols) <= i: + cols.append([]) + cols[i].extend(col) + n = len(cols[0]) if cols else 0 + rows = [[cols[c][i] if c < len(cols) else None for c in range(len(names))] + for i in range(n)] + return names, rows + + +def scalar(cypher: str, **kw) -> float | int | str | None: + _, rows = query(cypher, **kw) + return rows[0][0] if rows and rows[0] else None + + +# -------------------------------------------------------------------------- +# the named queries +# -------------------------------------------------------------------------- + +#: Named queries, each a (description, cypher) pair. These are the same +#: questions the studio offers; `crystal run --cypher` prints the Cypher +#: so it can be copied, edited and re-run with `crystal query`. +#: Named queries. Each is (description, cypher, tally-by). +#: +#: `tally-by` is not decoration. This Cypher dialect has no GROUP BY, and an +#: aggregate cannot be combined with another return item -- `RETURN +#: f.fragment_type, count(c)` is rejected, and so is an aggregate inside +#: ORDER BY. So a "how many of each" question is expressed as a FLAT +#: projection and counted client-side. When tally-by is set, the CLI groups on +#: that column and prints counts; `--cypher` still shows exactly the query that +#: ran, so nothing is hidden. +LIBRARY: dict[str, tuple[str, str, str | None]] = { + "spacegroups": ( + "Structures per space group, most common first", + "MATCH (s:Structure)-[:IN_SPACE_GROUP]->(g:SpaceGroup) " + "RETURN g.hm_symbol", + "g.hm_symbol", + ), + "centrosymmetric": ( + "How many structures are in a centrosymmetric space group", + "MATCH (s:Structure)-[:IN_SPACE_GROUP]->(g:SpaceGroup) " + "WHERE g.is_centrosymmetric = true RETURN count(s)", + None, + ), + "contacts": ( + "The contact populations, counted by kind and whether H was located", + "MATCH ()-[r:CONTACT]->() RETURN r.kind, r.h_inferred", + "r.kind+r.h_inferred", + ), + "hbond-geometry": ( + "Mean H...A over located-hydrogen bonds only", + "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' " + "AND r.h_inferred = false RETURN avg(r.length)", + None, + ), + "dimensionality": ( + "Periodic-net dimensionality of the hydrogen-bond net", + "MATCH (s:Structure) WHERE s.has_net_dim = true RETURN s.net_dim", + "s.net_dim", + ), + "dimer": ( + "R2,2(8) carboxylic-acid dimers (involution-generated)", + "MATCH (f1:Fragment)<-[:IN_FRAGMENT]-(a1:Atom)" + "-[h:CONTACT]->(a2:Atom)-[:IN_FRAGMENT]->(f2:Fragment) " + "WHERE f1.fragment_type = 'carboxylic_acid' " + "AND f2.fragment_type = 'carboxylic_acid' " + "AND h.kind = 'hbond' AND h.h_inferred = false " + "AND h.is_involution = true " + "AND a1.element = 'O' AND a2.element = 'O' " + "RETURN a1.cod_id, a1.label, a2.label, h.symop_triplet, h.length, h.angle", + None, + ), + "catemer": ( + "The same acid-acid bonds generated by a screw or glide: C(4) chains", + "MATCH (f1:Fragment)<-[:IN_FRAGMENT]-(a1:Atom)" + "-[h:CONTACT]->(a2:Atom)-[:IN_FRAGMENT]->(f2:Fragment) " + "WHERE f1.fragment_type = 'carboxylic_acid' " + "AND f2.fragment_type = 'carboxylic_acid' " + "AND h.kind = 'hbond' AND h.h_inferred = false " + "AND h.is_involution = false " + "AND a1.element = 'O' AND a2.element = 'O' " + "RETURN a1.cod_id, h.symop_triplet, h.length, h.angle", + None, + ), + "synthon-homo": ( + "Acid-to-acid O-H...O bonds (the homosynthon)", + "MATCH (f1:Fragment)<-[:IN_FRAGMENT]-(a1:Atom)" + "-[h:CONTACT]->(a2:Atom)-[:IN_FRAGMENT]->(f2:Fragment) " + "WHERE f1.fragment_type = 'carboxylic_acid' " + "AND f2.fragment_type = 'carboxylic_acid' " + "AND h.kind = 'hbond' AND h.h_inferred = false " + "AND a1.element = 'O' AND a2.element = 'O' RETURN count(h)", + None, + ), + "synthon-hetero": ( + "Acid-to-pyridine O-H...N bonds (the heterosynthon)", + "MATCH (f1:Fragment)<-[:IN_FRAGMENT]-(a1:Atom)" + "-[h:CONTACT]->(a2:Atom)-[:IN_FRAGMENT]->(f2:Fragment) " + "WHERE f1.fragment_type = 'carboxylic_acid' " + "AND f2.fragment_type = 'pyridine_nitrogen' " + "AND h.kind = 'hbond' AND h.h_inferred = false " + "AND a1.element = 'O' AND a2.element = 'N' RETURN count(h)", + None, + ), + "fragments": ( + "Molecules carrying each perceived functional group", + "MATCH (c:Component)-[:HAS_FRAGMENT]->(f:Fragment) " + "RETURN f.fragment_type", + "f.fragment_type", + ), + "molecule-atoms": ( + "Every atomic site belonging to one molecule, across all its structures", + "MATCH (c:Component)<-[:IN_COMPONENT]-(a:Atom) " + "WHERE c.inchikey = 'WPYMKLBDIGXBTP-UHFFFAOYSA-N' " + "RETURN a.cod_id, a.label, a.element", + None, + ), + "packing": ( + "Every stored contact of one structure, with its symmetry operation", + "MATCH (a:Atom)-[r:CONTACT]->(b:Atom) WHERE a.cod_id = 2229029 " + "RETURN a.label, b.label, r.kind, r.symop, r.length, r.angle", + None, + ), + "solvates": ( + "Which solvents appear, and in how many structures", + "MATCH (s:Structure)-[:CONTAINS_COMPONENT]->(c:Component) " + "WHERE c.is_solvent = true RETURN c.formula", + "c.formula", + ), + "recurrence": ( + "Molecules appearing in the most structures", + "MATCH (s:Structure)-[:CONTAINS_COMPONENT]->(c:Component) " + "WHERE c.is_solvent = false AND c.has_inchikey = true " + "RETURN c.inchikey", + "c.inchikey", + ), +} + + +# -------------------------------------------------------------------------- +# output +# -------------------------------------------------------------------------- + + +#: Columns whose integers are IDENTIFIERS, not quantities. A COD id printed as +#: "2,206,537" cannot be pasted into a search box or a script, so digit +#: grouping is applied to counts and never to these. +_ID_COLUMN = ("cod_id", "id", "number", "year", "atomic_number") + + +def is_id_column(name: str) -> bool: + leaf = name.rsplit(".", 1)[-1].strip("`") + return leaf in _ID_COLUMN + + +def fmt(v: object, *, group: bool = True) -> str: + if v is None: + return "-" + if isinstance(v, bool): + return "true" if v else "false" + if isinstance(v, float): + return f"{v:.4f}".rstrip("0").rstrip(".") + if isinstance(v, int): + return f"{v:,}" if group else str(v) + return str(v) + + +def print_table(names: list[str], rows: list[list], limit: int | None = None) -> None: + if not names: + print("(no columns)") + return + shown = rows[:limit] if limit else rows + group = [not is_id_column(n) for n in names] + cells = [[fmt(v, group=group[i] if i < len(group) else True) + for i, v in enumerate(r)] for r in shown] + widths = [len(n) for n in names] + for r in cells: + for i, c in enumerate(r): + if i < len(widths): + widths[i] = max(widths[i], len(c)) + print(" ".join(n.ljust(widths[i]) for i, n in enumerate(names))) + print(" ".join("-" * w for w in widths)) + for r in cells: + print(" ".join(c.ljust(widths[i]) for i, c in enumerate(r))) + if limit and len(rows) > limit: + print(f"... {len(rows) - limit:,} more rows (use --limit 0 for all)") + print(f"\n{len(rows):,} row(s)") + + +def print_csv(names: list[str], rows: list[list]) -> None: + w = csv.writer(sys.stdout) + w.writerow(names) + for r in rows: + w.writerow(["" if v is None else v for v in r]) + + +def tally_rows(names: list[str], rows: list[list], + column: str) -> tuple[list[str], list[list]]: + """Count rows by one column, most frequent first. + + The dialect has no GROUP BY, so grouping happens here. Doing it in Python + over a flat projection is exact -- it is the same rows the server would + have grouped -- and it keeps the shipped Cypher something you can paste + into any other client. + """ + from collections import Counter + + wanted = column.split("+") + idx = [names.index(c) for c in wanted if c in names] or [0] + counts: Counter[tuple] = Counter( + tuple("-" if r[i] is None else r[i] for i in idx) for r in rows) + ordered = sorted(counts.items(), key=lambda kv: (-kv[1], str(kv[0]))) + return ([names[i] for i in idx] + ["count"], + [list(k) + [v] for k, v in ordered]) + + +def emit(names, rows, args) -> None: + if args.json: + json.dump([dict(zip(names, r)) for r in rows], sys.stdout, indent=2, default=str) + print() + elif args.csv: + print_csv(names, rows) + else: + print_table(names, rows, None if args.limit == 0 else args.limit) + + +# -------------------------------------------------------------------------- +# commands +# -------------------------------------------------------------------------- + + +def cmd_facts(args) -> int: + g, h = args.graph, args.host + q = lambda c: scalar(c, graph=g, host=h) # noqa: E731 + print(f"graph: {g} host: {h}\n") + rows = [ + ("structures", q("MATCH (s:Structure) RETURN count(s)")), + ("atomic sites", q("MATCH (n:Atom) RETURN count(n)")), + ("components", q("MATCH (n:Component) RETURN count(n)")), + ("covalent bonds", q("MATCH ()-[r:BONDED_TO]->() RETURN count(r)")), + ("space groups", q("MATCH (n:SpaceGroup) RETURN count(n)")), + ] + for k, v in rows: + print(f" {k:<26}{fmt(v):>14}") + + print("\ncontact layer (three populations, never blended):") + pops = [ + ("strong H-bond (H located)", + "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' AND r.h_inferred = false RETURN count(r)"), + ("inferred (no H refined)", + "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' AND r.h_inferred = true RETURN count(r)"), + ("weak C-H...A", + "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond_weak' RETURN count(r)"), + ("halogen bond", + "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'halogen' RETURN count(r)"), + ("close_contact (excluded)", + "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'close_contact' RETURN count(r)"), + ] + for k, c in pops: + print(f" {k:<26}{fmt(q(c)):>14}") + + ha = q("MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' AND r.h_inferred = false RETURN avg(r.length)") + an = q("MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' AND r.h_inferred = false RETURN avg(r.angle)") + print(f"\n mean H...A {float(ha):.3f} A at {float(an):.2f} deg " + f"(strong, located-H population only)") + + total = q("MATCH (s:Structure) RETURN count(s)") + centro = q("MATCH (s:Structure)-[:IN_SPACE_GROUP]->(g:SpaceGroup) " + "WHERE g.is_centrosymmetric = true RETURN count(s)") + print(f" centrosymmetric {100 * float(centro) / float(total):.1f}% of structures") + + print("\nperiodic-net dimensionality (h_inferred excluded):") + lab = {0: "0D finite motif", 1: "1D chain", 2: "2D sheet", 3: "3D framework"} + scored = sum(float(q(f"MATCH (s:Structure) WHERE s.net_dim = {k} RETURN count(s)")) + for k in range(4)) + for k in range(4): + n = float(q(f"MATCH (s:Structure) WHERE s.net_dim = {k} RETURN count(s)")) + print(f" {lab[k]:<26}{fmt(int(n)):>10} {100 * n / max(scored, 1):>5.1f}%") + print(f" {'scored':<26}{fmt(int(scored)):>10}") + return 0 + + +def cmd_list(args) -> int: + width = max(len(k) for k in LIBRARY) + for name, (desc, _c, _t) in LIBRARY.items(): + print(f" {name.ljust(width)} {desc}") + print(f"\n{len(LIBRARY)} named queries. " + f"`crystal run --cypher` prints the Cypher.") + return 0 + + +def cmd_run(args) -> int: + if args.name not in LIBRARY: + near = [k for k in LIBRARY if args.name in k] + print(f"no query named '{args.name}'." + + (f" Did you mean: {', '.join(near)}?" if near else + " Run `crystal list`."), file=sys.stderr) + return 2 + desc, cypher, tally = LIBRARY[args.name] + if args.cypher: + print(cypher) + return 0 + if not (args.json or args.csv): + print(f"{desc}\n") + names, rows = query(cypher, graph=args.graph, host=args.host) + if tally: + names, rows = tally_rows(names, rows, tally) + emit(names, rows, args) + return 0 + + +def cmd_query(args) -> int: + names, rows = query(args.cypher, graph=args.graph, host=args.host) + emit(names, rows, args) + return 0 + + +def cmd_at(args) -> int: + """Run a query against the corpus as it stood at the end of one year.""" + _, ledger = query( + "MATCH (s:Snapshot) RETURN s.year, s.commit_hash ORDER BY s.year", + graph=LEDGER_GRAPH, host=args.host) + hashes = {int(y): h for y, h in ledger} + if args.year not in hashes: + print(f"no commit for {args.year}. Available: " + f"{min(hashes)}-{max(hashes)}", file=sys.stderr) + return 2 + commit = hashes[args.year] + # A commit must be resident before it can be read; this is idempotent. + query(f"LOAD COMMIT '{commit}'", graph=VERSIONED_GRAPH, host=args.host) + if not (args.json or args.csv): + print(f"cod@{args.year} commit {commit}\n") + names, rows = query(args.cypher, graph=VERSIONED_GRAPH, + host=args.host, commit=commit) + emit(names, rows, args) + return 0 + + +def cmd_commits(args) -> int: + names, rows = query( + "MATCH (s:Snapshot) RETURN s.year, s.tag, s.structures_total, " + "s.structures_added, s.commit_hash ORDER BY s.year", + graph=LEDGER_GRAPH, host=args.host) + emit(names, rows, args) + return 0 + + +#: The graphs this project builds. A server may carry others -- scratch graphs, +#: earlier builds -- and a new user has no way to tell which are the deliverable. +PROJECT_GRAPHS = ("cod_slice_v2", "cod_contacts_v2", "cod_versioned", "cod_versions") + + +def cmd_graphs(args) -> int: + names, rows = query("LIST AVAILABLE GRAPHS", graph="default", host=args.host) + if not (args.json or args.csv) and names: + names = list(names) + ["project"] + rows = [list(r) + ["yes" if r[0] in PROJECT_GRAPHS else "-"] + for r in rows] + rows.sort(key=lambda r: (r[-1] != "yes", str(r[0]))) + emit(names, rows, args) + if not (args.json or args.csv): + missing = [g for g in PROJECT_GRAPHS + if g not in {str(r[0]) for r in rows}] + if missing: + print(f"\nmissing: {', '.join(missing)} -- run ./run.sh --ingest") + return 0 + + +def main(argv: list[str] | None = None) -> int: + def add_common(parser: argparse.ArgumentParser, *, suppress: bool) -> None: + """Common flags, accepted either before OR after the subcommand. + + `crystal run dimer --limit 5` is the obvious thing to type, so it has + to work as well as `crystal --limit 5 run dimer`. The subcommand copies + default to SUPPRESS, so they only appear in the namespace when actually + given and therefore never clobber a value set before the subcommand. + """ + d = (lambda v: argparse.SUPPRESS) if suppress else (lambda v: v) + parser.add_argument("--host", default=d(DEFAULT_HOST), + help=f"TuringDB REST endpoint (default {DEFAULT_HOST})") + parser.add_argument("--graph", default=d(DEFAULT_GRAPH), + help=f"graph to query (default {DEFAULT_GRAPH})") + parser.add_argument("--csv", action="store_true", + default=d(False), help="CSV to stdout") + parser.add_argument("--json", action="store_true", + default=d(False), help="JSON to stdout") + parser.add_argument("--limit", type=int, default=d(40), + help="rows to print, 0 for all (default 40)") + + ap = argparse.ArgumentParser( + prog="crystal", description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + add_common(ap, suppress=False) + + common = argparse.ArgumentParser(add_help=False) + add_common(common, suppress=True) + + sub = ap.add_subparsers(dest="cmd", required=True) + + sub.add_parser("facts", parents=[common], + help="corpus summary").set_defaults(fn=cmd_facts) + sub.add_parser("list", parents=[common], + help="the named queries").set_defaults(fn=cmd_list) + sub.add_parser("commits", parents=[common], + help="the version ledger").set_defaults(fn=cmd_commits) + sub.add_parser("graphs", parents=[common], + help="graphs on the server").set_defaults(fn=cmd_graphs) + + p = sub.add_parser("run", parents=[common], help="run a named query") + p.add_argument("name") + p.add_argument("--cypher", action="store_true", help="print the Cypher, do not run") + p.set_defaults(fn=cmd_run) + + p = sub.add_parser("query", parents=[common], help="run your own Cypher") + p.add_argument("cypher") + p.set_defaults(fn=cmd_query) + + p = sub.add_parser("at", parents=[common], + help="run Cypher against the corpus as of a year") + p.add_argument("year", type=int) + p.add_argument("cypher") + p.set_defaults(fn=cmd_at) + + args = ap.parse_args(argv) + try: + return args.fn(args) + except QueryError as exc: + print(f"error: {exc}", file=sys.stderr) + return 1 + except BrokenPipeError: + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/ingest/build_graph.py b/ingest/build_graph.py index 8d56441..c17d5a7 100644 --- a/ingest/build_graph.py +++ b/ingest/build_graph.py @@ -24,12 +24,11 @@ import argparse import csv import json -import os +import math import signal import sys import time from collections import Counter -from dataclasses import dataclass, field from pathlib import Path from typing import Any, Iterable, Sequence @@ -40,7 +39,9 @@ from .contacts import detect_contacts from .download import cod_cif_relpath from .fragments import FRAGMENTS_BY_NAME -from .perception import perceive_component +from .netdim import is_involution as _involution_uncached +from .netdim import net_dimensionality +from .perception import perceive_component, perception_timeouts from .periodic import PeriodicNeighbours, asymmetric_positions from .solvents import assign_roles, solvent_name from .structure import ( @@ -52,6 +53,37 @@ read_structure, ) +#: Elements allowed to be an INFERRED hydrogen-bond donor, i.e. one where no +#: hydrogen was refined. Halogens are deliberately absent: see the note at the +#: contact-emission site. +INFERRED_DONOR_ELEMENTS = frozenset({"N", "O", "S"}) + +#: The symop alphabet is tiny (a few hundred distinct triplets over the whole +#: corpus) while the contact count is in the millions, so both of these are +#: memoised rather than recomputed per edge. +_INVOLUTION_CACHE: dict[str, bool] = {} +_CENTRO_CACHE: dict[int, bool] = {} + + +def _involution(triplet: str) -> bool: + got = _INVOLUTION_CACHE.get(triplet) + if got is None: + got = _involution_uncached(triplet) + _INVOLUTION_CACHE[triplet] = got + return got + + +def _centrosymmetric(number: int) -> bool: + got = _CENTRO_CACHE.get(number) + if got is None: + import gemmi as _g + + sg = _g.find_spacegroup_by_number(number) + got = bool(sg.is_centrosymmetric()) if sg is not None else False + _CENTRO_CACHE[number] = got + return got + + MAX_SITES = 400 # skip pathologically large asymmetric units #: Wall-clock budget for one structure. COD contains a long tail of entries that #: are individually pathological -- unusual symmetry, extreme disorder, or a @@ -146,6 +178,7 @@ def process_structure(args: tuple[str, int]) -> dict[str, Any] | None: "n_atoms": p.n_atoms, "n_heavy": p.n_heavy, "inchikey": p.inchikey, "fallback_key": p.fallback_key, "charge": p.charge, "fragments": p.fragments, + "fragment_sites": p.fragment_sites, "is_polymeric": False, "error": p.error, }) @@ -290,6 +323,53 @@ def edge(label: str, src: int, tgt: int, props: dict[str, Any] | None = None) -> if value is None: missing[field_name] += 1 + # ---- periodic-net dimensionality -------------------------------- + # Computed at emit rather than in the parse worker: it needs only what + # the parse cache already holds, so `--from-cache` can recompute it + # without re-parsing the corpus. + # + # THE COVALENT BONDS MUST BE IN THE GRAPH, not only the contacts. A + # hydrogen bond runs donor-heavy-atom -> acceptor-heavy-atom, and those + # two sites are joined in the crystal *through the molecules they belong + # to* -- an edge that is not a contact. Built from contacts alone the + # quotient graph is bipartite between donors and acceptors with no path + # longer than one hop: no cycles, rank 0, and almost every structure + # reports as a 0D finite motif, including textbook catemers. + # + # Including the intramolecular bonds contracts each molecule to a point, + # which is the standard formulation. It is safe: a ring INSIDE a finite + # molecule closes with the zero lattice vector and contributes nothing + # to the rank. + net_bonds = [(b["a"], b["b"], b["symop"]) for b in res["bonds"]] + uid_site = {a["uid"]: i for i, a in enumerate(res["atoms"])} + sg_no = res["spacegroup_number"] + + def _dim(kinds: frozenset[str], _res: dict[str, Any] = res, + _sg: int = sg_no, _bonds: list = net_bonds, + _sites: dict[str, int] = uid_site) -> int: + # The loop variables are bound as defaults rather than captured: + # the closure is called within the same iteration, but binding is + # correct by construction and keeps the intent explicit. + edges = [ + (c["a"], c["b"], c["symop"]) + for c in _res["contacts"] + if c["kind"] in kinds and not c["h_inferred"] + ] + if not edges: + return -1 # no qualifying contact; NOT a 0D claim + try: + d = net_dimensionality(_sg, _bonds + edges, _sites)[0] + except Exception: # noqa: BLE001 + return -2 # declined (over the expansion ceiling) + return d + + net_dim = _dim(frozenset({"hbond"})) + net_dim_weak = _dim(frozenset({"hbond", "hbond_weak"})) + # A coordination polymer's COVALENT net already carries lattice + # translations, so for those structures the number describes the + # framework, not the hydrogen-bond net. Flagged rather than mixed in. + net_polymeric = any(c["is_polymeric"] for c in res["components"]) + s_props: dict[str, Any] = { "cod_id": cod_id, "a": a, "b": b, "c": c, "alpha": al, "beta": be, "gamma": ga, @@ -307,6 +387,15 @@ def edge(label: str, src: int, tgt: int, props: dict[str, Any] | None = None) -> "has_z_prime": zprime is not None, "n_components": len(res["components"]), "n_contacts": len(res["contacts"]), + # -1 means "no qualifying contact", which is NOT the same claim as + # 0D. Kept distinct so a query can exclude it rather than count it + # as an isolated motif. + "net_dim": net_dim, + "net_dim_weak": net_dim_weak, + "has_net_dim": net_dim >= 0, + # true where the covalent net itself is periodic, so net_dim + # describes the framework rather than the hydrogen-bond net + "net_dim_is_framework": net_polymeric, # honest flag: this structure's contact search hit its time budget, # so its contact list is incomplete "contacts_truncated": bool(res.get("stats", {}).get("timed_out")), @@ -321,6 +410,10 @@ def edge(label: str, src: int, tgt: int, props: dict[str, Any] | None = None) -> "hm_symbol": res["spacegroup_hm"], "number": res["spacegroup_number"], "crystal_system": res["crystal_system"], + # Asked of gemmi rather than hand-listed: the centrosymmetric + # groups have a long tail, so any short hand-picked list + # understates the true share by several points. + "is_centrosymmetric": _centrosymmetric(res["spacegroup_number"]), }) edge("IN_SPACE_GROUP", sid, gid) @@ -355,23 +448,50 @@ def edge(label: str, src: int, tgt: int, props: dict[str, Any] | None = None) -> }) report["bonds"] += len(res["bonds"]) + element_of = {a["uid"]: a["element"] for a in res["atoms"]} for ct in res["contacts"]: x, y = atom_ids.get(ct["a"]), atom_ids.get(ct["b"]) if x is None or y is None: continue + kind = ct["kind"] + + # The h_inferred fallback fires per COMPONENT, so in a structure + # that is otherwise fully H-refined every heavy atom of a + # hydrogen-free fragment -- a perchlorate, a chloride, a fluorinated + # ring -- would qualify as a "donor" at D...A < 3.5 A with no + # angular test. Organic fluorine and covalently bound chlorine do + # not donate a proton. Only N, O and S are kept as inferred donors; + # the rest are still stored, but as close_contact, which no + # hydrogen-bond statistic touches. + if kind == "hbond" and ct["h_inferred"]: + donor_el = element_of.get(ct["donor"] or ct["a"], "") + if donor_el not in INFERRED_DONOR_ELEMENTS: + kind = "close_contact" + edge("CONTACT", x, y, { - "kind": ct["kind"], "length": ct["length"], + "kind": kind, "length": ct["length"], "angle": ct["angle"], "symop": ct["symop"], "symop_triplet": ct["triplet"], "h_inferred": ct["h_inferred"], "has_angle": ct["angle"] is not None, + # Is the generating operation its own inverse? This is what + # separates a RING motif from a CATEMER: an inversion, mirror or + # two-fold relates exactly two molecules, while a 2(1) screw or + # a glide has infinite order and generates a chain. Both have a + # two-fold rotation part, so nothing short of the exact + # composition distinguishes them -- and without it an R2,2(8) + # dimer query returns C(4) catemers. + "is_involution": _involution(ct["triplet"]), "cod_id": cod_id, }) - report[f"contact_{ct['kind']}"] += 1 - if ct["h_inferred"]: + report[f"contact_{kind}"] += 1 + if ct["h_inferred"] and kind == "hbond": report["contact_hbond_inferred"] += 1 report["contacts"] += len(res["contacts"]) # ---- Components + Fragments ------------------------------------- + # index within this structure -> the (corpus-wide, deduplicated) + # Component node id, so the atoms can be joined to their molecule below + comp_cid: dict[int, int] = {} for comp in res["components"]: key = comp["inchikey"] or comp["fallback_key"] or f"poly_{comp['formula']}" props = { @@ -398,6 +518,7 @@ def edge(label: str, src: int, tgt: int, props: dict[str, Any] | None = None) -> if comp["error"]: perception_errors[comp["error"].split(":")[0]] += 1 cid = reg.node("Component", key, props) + comp_cid[comp["index"]] = cid report["components"] += 1 edge("CONTAINS_COMPONENT", sid, cid, { "role": comp["role"], "n_sites": len(comp["sites"]), @@ -410,7 +531,42 @@ def edge(label: str, src: int, tgt: int, props: dict[str, Any] | None = None) -> }) if not _edge_exists(edges, cid, fid, "HAS_FRAGMENT"): edge("HAS_FRAGMENT", cid, fid) - report["fragment_links"] += 1 + report["fragment_links"] += 1 + + # Atom -> Fragment. HAS_FRAGMENT says a MOLECULE contains a + # group somewhere; this says WHICH ATOMS are in it. Without it + # a motif query can only ask "an oxygen of a molecule that has + # a carboxylic acid", which in a polyfunctional molecule also + # matches its alcohols, esters and nitro oxygens. + for site in comp.get("fragment_sites", {}).get(frag, ()): + at_uid = res["atoms"][site]["uid"] if site < len(res["atoms"]) else None + aid2 = atom_ids.get(at_uid) if at_uid else None + if aid2 is None: + continue + if not _edge_exists(edges, aid2, fid, "IN_FRAGMENT"): + edge("IN_FRAGMENT", aid2, fid) + report["in_fragment"] += 1 + + # ---- Atom -> Component ------------------------------------------ + # Without this edge there is no path from a CONTACT between two atoms to + # the MOLECULES those atoms belong to, which is precisely what a synthon + # question needs: "two carboxylic-acid components joined by a reciprocal + # pair of hydrogen bonds" is unanswerable when Atom and Component share + # no edge. The parser already knows the assignment -- it is what + # find_components returned -- so this is exact, not inferred, and free. + # + # Note the asymmetry it creates, which is the point: a Component node is + # deduplicated corpus-wide by identity, so one Component collects + # IN_COMPONENT edges from every structure it appears in. That is the + # "same molecule in forty structures is one node with forty edges" claim + # made traversable from the atomic level. + for at in res["atoms"]: + aid = atom_ids.get(at["uid"]) + cid2 = comp_cid.get(at.get("comp", -1)) + if aid is None or cid2 is None: + continue + edge("IN_COMPONENT", aid, cid2) + report["in_component"] += 1 # ---- Publication / Journal / Author ------------------------------ doi = (meta.get("doi") or "").strip() @@ -468,12 +624,47 @@ def edge(label: str, src: int, tgt: int, props: dict[str, Any] | None = None) -> n_supersedes += 1 report["supersedes"] = n_supersedes + # A handful of COD entries carry a coordinate or cell parameter that comes + # through as NaN. Python's json.dumps emits a bare `NaN`, which is NOT valid + # JSON, and TuringDB's loader rejects the whole 3.5 GB file on the first one + # -- `[json.exception.parse_error.101] ... last read: '"fract_x":N'` -- so a + # single bad atom costs the entire corpus. + # + # Dropping the record is not an option: LOAD JSONL requires node ids and + # relationship ids to be dense from 0 with no gaps, so removing one + # renumbers everything after it. The non-finite value is nulled instead, + # which simply leaves that property absent on that entity, and the count is + # reported rather than swallowed. + n_sanitised = 0 + + def _strip_nonfinite(obj: dict) -> dict: + props = obj.get("properties") + if not isinstance(props, dict): + return obj + for k, v in list(props.items()): + if isinstance(v, float) and not math.isfinite(v): + props[k] = None + return obj + + def _dump(obj: dict) -> str: + nonlocal n_sanitised + try: + return json.dumps(obj, separators=(",", ":"), allow_nan=False) + except ValueError: + n_sanitised += 1 + return json.dumps(_strip_nonfinite(obj), separators=(",", ":"), + allow_nan=False) + out_path.parent.mkdir(parents=True, exist_ok=True) with out_path.open("w") as fh: for node in reg.nodes: - fh.write(json.dumps(node, separators=(",", ":")) + "\n") + fh.write(_dump(node) + "\n") for e in edges: - fh.write(json.dumps(e, separators=(",", ":")) + "\n") + fh.write(_dump(e) + "\n") + report["nonfinite_sanitised"] = n_sanitised + if n_sanitised: + print(f" NON-FINITE : {n_sanitised:,} entities carried a NaN/inf " + f"property; those properties were nulled (JSON has no NaN)") summary = { "nodes": len(reg.nodes), @@ -488,13 +679,27 @@ def edge(label: str, src: int, tgt: int, props: dict[str, Any] | None = None) -> return summary +#: Edge types that must appear at most once between a given pair. Component, +#: Fragment, Journal and Author nodes are all deduplicated corpus-wide, so the +#: same pair is reached again every time another structure containing that +#: molecule, group, journal or author is emitted. +_SEEN_EDGES: set[tuple[str, int, int]] = set() + + def _edge_exists(edges: list[dict[str, Any]], src: int, tgt: int, label: str) -> bool: - # Cheap guard for the few edge types that are naturally deduplicated. - # Linear scan is fine because it is only consulted for Fragment/Journal. - s, t = str(src), str(tgt) - for e in reversed(edges[-400:]): - if e["label"] == label and e["start"]["id"] == s and e["end"]["id"] == t: - return True + """True if this (label, src, tgt) edge has already been emitted. + + Exact, via a set. A bounded scan of recent edges is not sufficient: a + Component node is shared by every structure containing that molecule, and + those structures are scattered through the corpus, so the second visit is + usually thousands of edges later. Under-deduplicating here multiplies any + query that joins Fragment on both ends of a contact -- the pair count + becomes (matches on side 1 x matches on side 2) rather than one. + """ + key = (label, src, tgt) + if key in _SEEN_EDGES: + return True + _SEEN_EDGES.add(key) return False @@ -641,15 +846,31 @@ def main(argv: Sequence[str] | None = None) -> int: rep_path.write_text(json.dumps(summary, indent=2) + "\n") c = summary["counts"] - print(f"\n=== ingest report ===") + print("\n=== ingest report ===") print(f" structures : {c.get('structures', 0):,} " f"(skipped {c.get('skipped', 0):,})") print(f" atoms : {c.get('atoms', 0):,}") print(f" bonds : {c.get('bonds', 0):,}") + print(f" Atom->Component : {c.get('in_component', 0):,}") + print(f" Atom->Fragment : {c.get('in_fragment', 0):,}") + # An empty Fragment layer is silent: perception failures do not raise, and + # the component and InChIKey counts are unaffected by them. Assert on it + # rather than expecting a reader to notice a missing line. + if c.get("components", 0) and not c.get("fragment_links", 0): + print(" !! FRAGMENT LAYER IS EMPTY -- perception is broken; every " + "synthon, coformer and solid-form query will return nothing") + _to = perception_timeouts() + if _to: + print(f" RDKit timeouts : {_to:,} components abandoned to a killed " + f"worker (see ingest/rdkit_guard.py)") print(f" contacts : {c.get('contacts', 0):,} " f"(hbond {c.get('contact_hbond', 0):,} of which " f"{c.get('contact_hbond_inferred', 0):,} h_inferred, " - f"halogen {c.get('contact_halogen', 0):,})") + f"weak C-H...A {c.get('contact_hbond_weak', 0):,}, " + f"halogen {c.get('contact_halogen', 0):,}, " + f"close_contact {c.get('contact_close_contact', 0):,} " + f"[inferred with a non-N/O/S donor, excluded from every " + f"hydrogen-bond statistic])") print(f" components : {c.get('components', 0):,} " f"(InChIKey {c.get('components_with_inchikey', 0):,}, " f"fallback {c.get('components_fallback', 0):,}, " diff --git a/ingest/chemistry.py b/ingest/chemistry.py index 74b2345..8c56d04 100644 --- a/ingest/chemistry.py +++ b/ingest/chemistry.py @@ -127,6 +127,32 @@ #: Elements admitted as hydrogen-bond donor or acceptor heavy atoms. HBOND_ELEMENTS: Final[frozenset[str]] = frozenset({"N", "O", "F", "S", "Cl"}) +#: Weak hydrogen bond: C-H...A, the interaction Taylor & Kennard established is +#: real and that Desiraju & Steiner made a structural tool. Admitting it matters +#: here for a specific reason: with only strong N/O-H...O/N bonds the stored +#: contact layer averages ~3.4 edges per structure and most packings do not +#: percolate, which is why the periodic-net dimensionality reads 0D for 88% of +#: the corpus. C-H...O is the cheapest correction -- same geometry, different +#: donor set, no ring perception needed. +#: +#: Two rules the detector enforces, both deliberate: +#: +#: * These edges are written with kind='hbond_weak', never 'hbond'. Every query +#: that filters kind='hbond' therefore returns exactly what it returned +#: before. The project rule that the strong and inferred populations are never +#: silently mixed now extends to three populations. +#: * There is NO h_inferred fallback for the weak class. A strong N...O at 3.5 A +#: with no located hydrogen is a defensible inference; a C...O at 3.5 A with no +#: located hydrogen is just two atoms near each other. Weak contacts are +#: emitted only where the hydrogen was actually refined. +#: +#: 2.90 A is the conventional H...A ceiling for a structurally meaningful +#: C-H...O (Steiner 2002 Table 2 puts the bulk of the population inside it); the +#: 120 deg angle floor is the same one used for the strong class. +WEAK_HBOND_DONORS: Final[frozenset[str]] = frozenset({"C"}) +WEAK_HBOND_MAX_H_ACCEPTOR: Final = 2.90 # angstrom +WEAK_HBOND_MIN_ANGLE: Final = 120.0 # degrees + #: Halogen bond: C-X...A where X is Cl, Br or I. The X...A distance must be #: below the sum of van der Waals radii, and the C-X...A angle above 150 deg #: (the sigma-hole is directed along the C-X axis). @@ -199,6 +225,10 @@ def is_halogen_donor(symbol: str) -> bool: return normalise_element(symbol) in HALOGEN_DONORS +def is_weak_hbond_donor(symbol: str) -> bool: + return normalise_element(symbol) in WEAK_HBOND_DONORS + + #: Reported in the ingest report so the exact conventions travel with the data. CRITERIA_PROVENANCE: Final[dict[str, str]] = { "covalent_radii": ( @@ -219,6 +249,15 @@ def is_halogen_donor(symbol: str) -> bool: "hbond_without_h": ( f"D...A < {HBOND_HEAVY_MAX} A, flagged h_inferred=true" ), + "hbond_weak": ( + f"C-H...A with a LOCATED H; H...A < {WEAK_HBOND_MAX_H_ACCEPTOR} A and " + f"C-H...A > {WEAK_HBOND_MIN_ANGLE} deg, A in {sorted(HBOND_ELEMENTS)}; " + "stored as kind='hbond_weak', never merged with the strong population, " + "and never inferred without a hydrogen. " + "Taylor & Kennard, J. Am. Chem. Soc., 1982, 104, 5063 " + "(doi:10.1021/ja00383a012); Desiraju & Steiner, The Weak Hydrogen Bond, " + "OUP 1999; Steiner, Angew. Chem. Int. Ed., 2002, 41, 48" + ), "halogen_bond": ( f"X in {sorted(HALOGEN_DONORS)}; X...A < sum of vdW radii and " f"C-X...A > {HALOGEN_MIN_ANGLE} deg" diff --git a/ingest/contacts.py b/ingest/contacts.py index a9a4437..fc28d16 100644 --- a/ingest/contacts.py +++ b/ingest/contacts.py @@ -30,7 +30,6 @@ from . import chemistry as ch from .periodic import PeriodicNeighbours from .structure import Atom, Bond, _op_for -from .symmetry import SymOp #: Bounds the per-component molecule walk. A component larger than this is an #: extended framework, which structure.find_components already flags. @@ -45,7 +44,7 @@ class Contact: b_index: int a_uid: str b_uid: str - kind: str # 'hbond' | 'halogen' | 'pi_stack' + kind: str # 'hbond' | 'hbond_weak' | 'halogen' | 'pi_stack' length: float # H...A, X...A, or centroid...centroid angle: float | None # D-H...A, C-X...A, or interplanar symop: str @@ -127,6 +126,12 @@ def detect_contacts( the criterion exists for, and correctly still covers an unrefined water in a structure whose organic component does have hydrogens. + Three populations come out of this function and they are kept apart by + `kind`: 'hbond' (strong, located H), 'hbond' with h_inferred=True (strong, + heavy-atom fallback), and 'hbond_weak' (C-H...A, located H only). A query + filtering kind='hbond' sees exactly what it saw before weak contacts + existed. + pi-stacking is not produced here: it needs ring perception, and the spec's centroid criterion alone admits badly offset stacks that most crystal engineers would not call pi-stacking. Shipping a half-defined criterion @@ -141,7 +146,8 @@ def detect_contacts( whole ingest on one entry. """ stats: dict[str, int] = { - "hbond": 0, "hbond_inferred": 0, "halogen": 0, "intramolecular_skipped": 0, + "hbond": 0, "hbond_inferred": 0, "hbond_weak": 0, "halogen": 0, + "intramolecular_skipped": 0, } out: list[Contact] = [] adj = build_adjacency(sg, bonds) @@ -231,10 +237,18 @@ def donor_molecule_lacks_hydrogen(site: int) -> bool: seen: set[tuple[int, int, str, str]] = set() search_r = max(ch.HBOND_HEAVY_MAX, ch.MAX_CONTACT_SEARCH_RADIUS) + # A weak-bond donor is a carbon that actually carries a refined hydrogen. + # Gating on `hydrogens` here rather than inside the loop is what keeps the + # cost down: in an entry with no H positions at all (much of Acta E) the + # carbon set collapses to empty and the weak pass costs nothing. donors = [ a for a in atoms if a.occupancy >= ch.MIN_BONDING_OCCUPANCY - and (ch.is_hbond_element(a.element) or ch.is_halogen_donor(a.element)) + and ( + ch.is_hbond_element(a.element) + or ch.is_halogen_donor(a.element) + or (ch.is_weak_hbond_donor(a.element) and a.site_index in hydrogens) + ) ] if not donors: return out, stats @@ -343,6 +357,42 @@ def donor_molecule_lacks_hydrogen(site: int) -> bool: stats["hbond"] += 1 stats["hbond_inferred"] += 1 + # ---- weak hydrogen bond, C-H...A ----------------------------- + # Same geometry as the strong class, a looser distance ceiling, and + # a separate `kind` so the two never merge in a query. No inferred + # variant: without a located H this is not a hydrogen bond, it is + # two atoms near each other (see chemistry.WEAK_HBOND_DONORS). + if ch.is_weak_hbond_donor(donor.element): + for h_site, h_op in hydrogens.get(donor.site_index, ()): + h_atom = by_index.get(h_site) + if h_atom is None: + continue + hf = h_op.apply_to_xyz( + [h_atom.fract_x, h_atom.fract_y, h_atom.fract_z] + ) + hp = st.cell.orthogonalize(gemmi.Fractional(*hf)) + h_pos = np.array([hp.x, hp.y, hp.z]) + d_ha = float(np.linalg.norm(a_pos - h_pos)) + if d_ha >= ch.WEAK_HBOND_MAX_H_ACCEPTOR: + continue + angle = _angle_deg(d_pos, h_pos, a_pos) + if angle <= ch.WEAK_HBOND_MIN_ANGLE: + continue + kw = (donor.site_index, acceptor.site_index, + symop.cif_code(), f"w{h_atom.uid}") + if kw in seen: + continue + seen.add(kw) + out.append(Contact( + a_index=donor.site_index, b_index=acceptor.site_index, + a_uid=donor.uid, b_uid=acceptor.uid, kind="hbond_weak", + length=round(d_ha, 4), angle=round(angle, 2), + symop=symop.cif_code(), + symop_triplet=symop.full_triplet(), h_inferred=False, + donor_uid=donor.uid, hydrogen_uid=h_atom.uid, + )) + stats["hbond_weak"] += 1 + # ---- halogen bond -------------------------------------------- if ch.is_halogen_donor(donor.element): cutoff = ch.vdw_sum(donor.element, acceptor.element) diff --git a/ingest/download.py b/ingest/download.py index c327995..698943c 100644 --- a/ingest/download.py +++ b/ingest/download.py @@ -26,7 +26,6 @@ import csv import hashlib import json -import os import re import shutil import subprocess @@ -37,7 +36,7 @@ from pathlib import Path from typing import Any, Iterable, Sequence -from .cod_journals import JOURNALS, JOURNALS_BY_KEY, all_variants, journal_key_for +from .cod_journals import JOURNALS, JOURNALS_BY_KEY, journal_key_for # -------------------------------------------------------------------------- # constants diff --git a/ingest/netdim.py b/ingest/netdim.py new file mode 100644 index 0000000..ff53579 --- /dev/null +++ b/ingest/netdim.py @@ -0,0 +1,218 @@ +"""Periodic-net dimensionality: the rank of the cycle lattice. + +Extracted from queries/dimensionality_net.py so that ingest/build_graph.py can +stamp the result onto every Structure node at build time. Having it as a node +property is what lets the UI ask "show me the 2D sheets" as a one-line MATCH +instead of shipping a Python script alongside the demo. + +Why the naive traversal is wrong +-------------------------------- +The contact edges we store form a *quotient* graph: nodes are asymmetric-unit +sites and each edge carries the symmetry operation that generated the +neighbour. A hydrogen-bonded chain running through the crystal is therefore NOT +a long path in that graph -- it is a short cycle, because after a few hops you +return to the same asymmetric-unit site one unit cell over. Counting reachable +nodes by depth saturates immediately and reports every structure as isolated. + +The correct method is the standard periodic-net calculation. Build the labelled +quotient graph of the full unit cell, where every edge carries the integer +lattice vector crossed when traversing it. Take a spanning tree; for every +non-tree edge compute the net lattice translation around the cycle it closes. +The **rank of the set of those translations** is the dimensionality: + + rank 0 -> 0D finite motif (dimer, ring, isolated cluster) + rank 1 -> 1D chain + rank 2 -> 2D sheet + rank 3 -> 3D framework + +This is precisely what the symop on every edge is for: it decides the +dimensionality of an infinite periodic network without ever materialising a +supercell. +""" + +from __future__ import annotations + +from collections import defaultdict + +import gemmi +import numpy as np + +LABELS = {0: "0D finite motif", 1: "1D chain", 2: "2D sheet", 3: "3D framework"} +CONSEQUENCE = { + 0: "isolated dimers/clusters - typically the most soluble", + 1: "chains - often needle morphology, anisotropic mechanical response", + 2: "layers - plate morphology, easy cleavage, good tabletting", + 3: "three-dimensional network - harder, less soluble, humidity stable", +} + + +def _compose(g: gemmi.Op, s: gemmi.Op) -> tuple[np.ndarray, np.ndarray]: + """Compose g after s, KEEPING the lattice part of the translation. + + gemmi's own `g * s` reduces the translation modulo one unit cell, which + silently discards exactly the information this calculation depends on: + + S = -x,-y+1,-z tran [0, 24, 0] + g[0] * S = -x,-y,-z tran [0, 0, 0] <- the +1 in y is gone + + With every lattice translation flattened to zero, every cycle closes with a + zero period and every structure looks 0D. So we do the composition in exact + integer arithmetic instead. gemmi scales both rot and tran by Op.DEN, so + with R = rot/DEN an integer rotation matrix: + + rot(g.s) = R_g @ rot_s (still DEN-scaled) + tran(g.s) = R_g @ tran_s + tran_g (still DEN-scaled, NOT reduced) + """ + den = gemmi.Op.DEN + rg = np.array(g.rot, dtype=np.int64) // den + rot = rg @ np.array(s.rot, dtype=np.int64) + tran = rg @ np.array(s.tran, dtype=np.int64) + np.array(g.tran, dtype=np.int64) + return rot, tran + + +def _decompose( + rot: np.ndarray, tran: np.ndarray, ops: list[gemmi.Op] +) -> tuple[int, tuple[int, int, int]] | None: + """Write (rot, tran) as ops[k] followed by an integer lattice translation.""" + den = gemmi.Op.DEN + for k, op in enumerate(ops): + if not np.array_equal(np.array(op.rot, dtype=np.int64), rot): + continue + d = tran - np.array(op.tran, dtype=np.int64) + if np.all(d % den == 0): + q = d // den + return k, (int(q[0]), int(q[1]), int(q[2])) + return None + + +def _op_from_code(code: str, ops: list[gemmi.Op]) -> gemmi.Op: + idx_s, tr_s = code.split("_") + base = ops[int(idx_s) - 1] + op = gemmi.Op(base.triplet()) + den = gemmi.Op.DEN + sh = [int(c) - 5 for c in tr_s] + op.tran = (op.tran[0] + sh[0] * den, + op.tran[1] + sh[1] * den, + op.tran[2] + sh[2] * den) + return op + + +def net_dimensionality( + spacegroup_number: int, + contacts: list[tuple[str, str, str]], + uid_to_site: dict[str, int], +) -> tuple[int, int, int]: + """Return (dimensionality, n_nodes_in_largest_component, n_components). + + `contacts` is [(a_uid, b_uid, symop_code)]; a_uid sits at the identity and + b_uid is generated by symop_code. + """ + sg = gemmi.find_spacegroup_by_number(spacegroup_number) + if sg is None or not contacts: + return 0, 0, 0 + ops = list(sg.operations()) + n_ops = len(ops) + # Hard ceiling on the orbit expansion. A hang here is disproportionately + # expensive -- it stalls a whole-corpus ingest with no error -- so an + # over-large structure is declined rather than attempted. + if len(contacts) * n_ops * n_ops > 4_000_000: + return -1, 0, 0 + + # Unit-cell nodes are (asymmetric site, symmetry operation). Applying every + # operation g to a stored edge gives the whole orbit of that contact. + adj: dict[tuple[int, int], list[tuple[tuple[int, int], np.ndarray]]] = defaultdict(list) + for a_uid, b_uid, code in contacts: + a = uid_to_site.get(a_uid) + b = uid_to_site.get(b_uid) + if a is None or b is None: + continue + try: + s_op = _op_from_code(code, ops) + except Exception: # noqa: BLE001 + continue + for g_idx, g in enumerate(ops): + c_rot, c_tran = _compose(g, s_op) + got = _decompose(c_rot, c_tran, ops) + if got is None: + continue + k, tau = got + v = np.array(tau, dtype=np.int64) + adj[(a, g_idx)].append(((b, k), v)) + adj[(b, k)].append(((a, g_idx), -v)) + + if not adj: + return 0, 0, 0 + + # spanning forest with accumulated lattice potentials; every non-tree edge + # closes a cycle whose net translation is a period of the net. + # + # Cycles are collected as a SET of integer triples, not a list. A dense net + # closes one cycle per non-tree edge -- hundreds of thousands of them on a + # well-connected structure once weak C-H...A contacts are admitted -- and + # feeding all of those to matrix_rank builds an enormous (n x 3) matrix + # whose rank can only ever be 0..3. Worse, it spends that time inside LAPACK + # with the GIL released, so the ingest's SIGALRM watchdog cannot fire and a + # single structure can hang the whole run at 100% CPU with no diagnostic. + # + # De-duplicating is exact -- rank depends on the SET of vectors, not on how + # many times each appears -- and collapses those hundreds of thousands to a + # handful of distinct primitive translations. + seen: dict[tuple[int, int], np.ndarray] = {} + cycle_set: set[tuple[int, int, int]] = set() + comp_sizes: list[int] = [] + for start in list(adj): + if start in seen: + continue + seen[start] = np.zeros(3, dtype=np.int64) + stack = [start] + size = 0 + while stack: + u = stack.pop() + size += 1 + for v, t in adj.get(u, ()): + if v not in seen: + seen[v] = seen[u] + t + stack.append(v) + else: + cyc = seen[u] + t - seen[v] + if cyc.any(): + cycle_set.add((int(cyc[0]), int(cyc[1]), int(cyc[2]))) + comp_sizes.append(size) + + if not cycle_set: + return 0, max(comp_sizes, default=0), len(comp_sizes) + rank = int(np.linalg.matrix_rank( + np.array(sorted(cycle_set), dtype=np.float64) + )) + return min(rank, 3), max(comp_sizes, default=0), len(comp_sizes) + + +def is_involution(triplet: str) -> bool: + """True if applying this symmetry operation twice is the identity EXACTLY. + + This is the discriminator between a ring motif and a catemer, and it is the + thing that decides whether an R2,2(8) acid dimer query is right or wrong. + + A centrosymmetric dimer is generated by an inversion, a mirror or a + two-fold: S o S is the identity with a ZERO lattice translation, the + operation has order 2, and it relates exactly two molecules. A catemer is + generated by a 2(1) screw, a glide or a pure translation: S o S is a + NON-ZERO lattice translation, the operation has infinite order, and it + generates a chain. + + Both have a two-fold rotation part, so testing the rotation alone does not + separate them -- the translation is the whole signal, which is why this + uses the exact integer composition from `_compose` rather than gemmi's + `g * s` (that reduces the translation modulo one cell and would call a + screw axis an involution). + """ + try: + op = gemmi.Op(triplet) + except Exception: # noqa: BLE001 + return False + rot, tran = _compose(op, op) + # _compose returns gemmi's DEN-scaled integers, so the identity rotation is + # DEN * I, not I. + den = gemmi.Op.DEN + identity = np.eye(3, dtype=np.int64) * den + return bool(np.array_equal(rot, identity) and not tran.any()) diff --git a/ingest/perception.py b/ingest/perception.py index 4f57b1d..1ef4314 100644 --- a/ingest/perception.py +++ b/ingest/perception.py @@ -18,11 +18,12 @@ from __future__ import annotations from collections import Counter -from dataclasses import dataclass +from dataclasses import dataclass, field import gemmi from .fragments import FragmentPattern, compiled_patterns +from .rdkit_guard import BondOrderPerceiver from .structure import Atom, Bond, hill_formula #: Above this atom count we do not attempt bond-order perception. RDKit's @@ -33,6 +34,23 @@ _PATTERNS: list[tuple[FragmentPattern, object]] | None = None +#: One killable RDKit worker per ingest process, created on first use. +_PERCEIVER: "BondOrderPerceiver | None" = None + + +def _perceiver() -> "BondOrderPerceiver": + global _PERCEIVER + if _PERCEIVER is None: + from .rdkit_guard import BondOrderPerceiver + + _PERCEIVER = BondOrderPerceiver() + return _PERCEIVER + + +def perception_timeouts() -> int: + """How many components were abandoned to a killed worker.""" + return _PERCEIVER.timeouts if _PERCEIVER is not None else 0 + def _has_metal(atoms: list[Atom], instances: list[tuple[int, gemmi.Op]]) -> bool: for si, _ in instances: @@ -60,6 +78,9 @@ class Perceived: n_heavy: int charge: int | None fragments: list[str] + #: fragment type -> the SITE INDICES it matched, so a query can ask which + #: atom belongs to the group rather than only which molecule contains one + fragment_sites: dict[str, list[int]] = field(default_factory=dict) error: str = "" @@ -72,7 +93,6 @@ def perceive_component( ) -> Perceived: """Build an RDKit molecule for one component and derive its identity.""" from rdkit import Chem, RDLogger - from rdkit.Chem import rdDetermineBonds from rdkit.Geometry import Point3D RDLogger.DisableLog("rdApp.*") @@ -146,17 +166,33 @@ def perceive_component( try: m = mol.GetMol() - # infer bond orders and formal charges from the 3D geometry; this is the - # step that fails on metal complexes - rdDetermineBonds.DetermineBondOrders(m, charge=0, embedChiral=False) - Chem.SanitizeMol(m) - inchi = Chem.MolToInchi(m) - if not inchi: - raise ValueError("empty InChI") - key = Chem.InchiToInchiKey(inchi) - charge = Chem.GetFormalCharge(m) - frags = _match_fragments(m) - return Perceived(key, fallback, formula, n_atoms, n_heavy, charge, frags) + # Infer bond orders and formal charges from the 3D geometry. This is the + # step that fails on metal complexes -- and, on a small number of + # ordinary organics, the step that NEVER RETURNS. It is a C++ call, so + # the ingest's own SIGALRM cannot interrupt it; it runs in a child + # process the parent can kill instead. See ingest/rdkit_guard.py. + status, payload = _perceiver().perceive(Chem.MolToMolBlock(m, kekulize=False)) + if status != "ok": + raise ValueError(str(payload)) + key = payload["inchikey"] + charge = payload["charge"] + # Fragment matching must run on the PERCEIVED molecule, not on `m`. + # Bond-order perception happens in the child process, so `m` still + # carries none and every bond-order-sensitive SMARTS would fail to + # match -- silently, since an unmatched pattern is indistinguishable + # from an absent group. The perceived structure comes back as SMILES: + # atom order is irrelevant to a SMARTS census, connectivity and bond + # order are not. + # The child matched the SMARTS on the perceived molecule, so these + # indices are in the caller's own atom order and map straight back to + # crystallographic sites through `instances`. + frags = payload.get("fragments") or [] + frag_sites = { + name: sorted({instances[i][0] for i in idxs if i < len(instances)}) + for name, idxs in (payload.get("fragment_atoms") or {}).items() + } + return Perceived(key, fallback, formula, n_atoms, n_heavy, charge, frags, + fragment_sites=frag_sites) except Exception as exc: # noqa: BLE001 # Still try fragment matching on the unsanitised connectivity: many # SMARTS here are purely topological and match without bond orders. diff --git a/ingest/periodic.py b/ingest/periodic.py index b08fa99..bf6c6b0 100644 --- a/ingest/periodic.py +++ b/ingest/periodic.py @@ -90,7 +90,6 @@ def __init__( # index means two genuinely distinct atoms are never merged, so the # tolerance only has to absorb refinement noise about the exact # special position. - occupied: set[tuple[int, int, int, int]] = set() inv_tol = 1.0 / SPECIAL_POSITION_TOLERANCE # bounding box of the asymmetric unit in fractional space, padded by the diff --git a/ingest/rdkit_guard.py b/ingest/rdkit_guard.py new file mode 100644 index 0000000..450ce0d --- /dev/null +++ b/ingest/rdkit_guard.py @@ -0,0 +1,162 @@ +"""A hard timeout around RDKit's bond-order perception. + +Why this exists +--------------- +`rdDetermineBonds.DetermineBondOrders` searches over bond-order assignments and +on some connectivity graphs it does not terminate. It is a C++ call, so it never +returns to the interpreter and **no Python-level watchdog can stop it**: the +ingest already arms a 20 s `SIGALRM` per structure and the signal simply sits +undelivered while the process spins at 100% CPU. + +An external `timeout` around the whole ingest does not solve it either: it +sends SIGTERM to a process that cannot act on one, so the child survives and +the wrapper is left supervising a process it cannot stop. + +The timeout therefore has to live in a process the parent can actually *kill*. + +Design +------ +One long-lived worker process, `fork`-ed rather than spawned: spawn re-imports +`__main__` in the child, which misbehaves when the caller is not an importable +module. Work goes over a `Pipe`, and the parent waits with a deadline: + +* answer inside the budget -> return it, keep the worker +* deadline passes -> SIGKILL the worker, respawn it lazily, report a + timeout for this component only + +The worker is reused across calls, so the cost is one pipe round trip per +component (measured well under a millisecond) rather than a process spawn. + +The molecule crosses the boundary as a mol block, which is the cheapest +representation RDKit will round-trip without needing the caller's RWMol. +""" + +from __future__ import annotations + +import multiprocessing as mp +import os +from typing import Any + +#: Budget for one component's bond-order perception. Generous next to the +#: sub-millisecond typical case, small next to "never returns". +PERCEPTION_TIMEOUT_S = 10.0 + + +def _worker(conn: Any) -> None: # pragma: no cover - runs in a child process + """Perceive mol blocks until the pipe closes. Killed, never asked to stop.""" + from rdkit import Chem, RDLogger + from rdkit.Chem import rdDetermineBonds + + RDLogger.DisableLog("rdApp.*") + + _compiled: list = [] + + def _patterns() -> list: + if not _compiled: + from .fragments import compiled_patterns + + _compiled.extend(compiled_patterns()) + return _compiled + + while True: + try: + molblock = conn.recv() + except EOFError: + return + try: + m = Chem.MolFromMolBlock(molblock, sanitize=False, removeHs=False) + if m is None: + conn.send(("error", "MolFromMolBlock returned None")) + continue + rdDetermineBonds.DetermineBondOrders(m, charge=0, embedChiral=False) + Chem.SanitizeMol(m) + # Fragment matching happens HERE, on the perceived molecule, while + # the atom order is still the one the caller built. The caller + # cannot do it: the only representation that crosses back is a + # SMILES string, whose atom order is canonical, not original -- so + # matched indices would be unmappable to crystallographic sites. + # Matching here is what makes atom-level fragment membership + # possible at all. + matches: dict[str, list[int]] = {} + for pattern, query in _patterns(): + try: + hits = m.GetSubstructMatches(query, uniquify=True) + except Exception: # noqa: BLE001 + continue + if hits: + matches[pattern.name] = sorted({i for h in hits for i in h}) + inchi = Chem.MolToInchi(m) + if not inchi: + conn.send(("error", "empty InChI")) + continue + conn.send(( + "ok", + { + "inchikey": Chem.InchiToInchiKey(inchi), + "charge": Chem.GetFormalCharge(m), + "fragments": sorted(matches), + "fragment_atoms": matches, + }, + )) + except Exception as exc: # noqa: BLE001 + conn.send(("error", type(exc).__name__ + ": " + str(exc)[:120])) + + +class BondOrderPerceiver: + """Perceive bond orders in a killable child process.""" + + def __init__(self, timeout_s: float = PERCEPTION_TIMEOUT_S) -> None: + self.timeout_s = timeout_s + self._proc: mp.process.BaseProcess | None = None + self._conn: Any = None + #: components abandoned because the child had to be killed + self.timeouts = 0 + + def _ensure(self) -> None: + if self._proc is not None and self._proc.is_alive(): + return + # fork, not spawn: the child must not re-import __main__, and the work + # is pure CPU holding no threads or open sockets. + ctx = mp.get_context("fork") + parent, child = ctx.Pipe() + self._proc = ctx.Process(target=_worker, args=(child,), daemon=True) + self._proc.start() + child.close() + self._conn = parent + + def _kill(self) -> None: + if self._proc is not None: + try: + if self._proc.pid: + os.kill(self._proc.pid, 9) + except (ProcessLookupError, PermissionError, AttributeError): + pass + self._proc.join(timeout=2) + self._proc = None + if self._conn is not None: + try: + self._conn.close() + except OSError: + pass + self._conn = None + + def perceive(self, molblock: str) -> tuple[str, Any]: + """Return ('ok', {...}) or ('error', reason). Never raises, never hangs.""" + try: + self._ensure() + self._conn.send(molblock) + if not self._conn.poll(self.timeout_s): + self._kill() + self.timeouts += 1 + return "error", ( + f"RDKitTimeout: bond-order perception exceeded " + f"{self.timeout_s:g}s and the worker was killed" + ) + return self._conn.recv() # type: ignore[no-any-return] + except (EOFError, BrokenPipeError, OSError) as exc: + # the child died mid-call; treat exactly like a failed perception + self._kill() + return "error", f"RDKitWorkerLost: {type(exc).__name__}: {exc}"[:160] + + def close(self) -> None: + self._kill() diff --git a/ingest/solvents.py b/ingest/solvents.py index 8a2bc2d..aa8251e 100644 --- a/ingest/solvents.py +++ b/ingest/solvents.py @@ -33,7 +33,6 @@ "C2 H6 O": "ethanol", "C2 O": "ethanol (H not located)", "C3 H8 O": "propanol / isopropanol", - "C3 O": "propanol (H not located)", "C2 H6 O S": "dimethyl sulfoxide", "C2 O S": "DMSO (H not located)", "C3 H7 N O": "dimethylformamide", @@ -41,7 +40,9 @@ "C2 H3 N": "acetonitrile", "C2 N": "acetonitrile (H not located)", "C3 H6 O": "acetone", - "C3 O": "acetone (H not located)", + # C3O with no located hydrogens is ambiguous between propan-1-ol, + # propan-2-ol and acetone; naming one of them would be a guess. + "C3 O": "propanol or acetone (H not located)", "C4 H8 O": "tetrahydrofuran", "C4 O": "THF (H not located)", "C4 H10 O": "diethyl ether", diff --git a/ingest/structure.py b/ingest/structure.py index 52d4b48..8368a8d 100644 --- a/ingest/structure.py +++ b/ingest/structure.py @@ -19,6 +19,8 @@ from __future__ import annotations +import math + import hashlib from collections import Counter, deque from dataclasses import dataclass, field @@ -29,7 +31,7 @@ from . import chemistry as ch from .periodic import PeriodicNeighbours, asymmetric_positions -from .symmetry import SymOp, resolve_symop, spacegroup_of +from .symmetry import SymOp, spacegroup_of #: Assembling a molecule stops here. Anything larger is an extended framework #: (coordination polymer, MOF, hydrogen-bonded network in the covalent graph), @@ -144,7 +146,7 @@ def _clean_float(value: float | None) -> float | None: f = float(value) except (TypeError, ValueError): return None - if f != f: # NaN + if not math.isfinite(f): return None return f diff --git a/pyproject.toml b/pyproject.toml index 4ca1da3..0d65153 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,3 +1,7 @@ +[build-system] +requires = ["hatchling"] +build-backend = "hatchling.build" + [project] name = "ccdc-demo" version = "0.1.0" @@ -12,10 +16,7 @@ dependencies = [ "rdkit>=2024.3.1", "numpy>=1.26", "scipy>=1.13", - # Metadata dump handling + ingest reporting. - "pandas>=2.2", - "pyarrow>=16.0", - # Downloader. + # Downloader: HTTP for the release manifest, progress for the long rsync. "httpx>=0.27", "tqdm>=4.66", # COD publishes a public read-only MySQL mirror; pure-python driver keeps @@ -24,31 +25,38 @@ dependencies = [ # TuringDB Python SDK. The server binary ships inside this wheel, so there # is no separate install step. "turingdb", - # Backend. - "fastapi>=0.111", - "uvicorn[standard]>=0.30", - "pydantic>=2.7", ] +[project.scripts] +crystal = "crystal:main" + [project.optional-dependencies] -bench = [ - "neo4j>=5.20", -] dev = [ - "pytest>=8.2", "mypy>=1.10", "ruff>=0.5", ] +[tool.hatch.build.targets.wheel] +packages = ["ingest", "queries"] +only-include = ["ingest", "queries", "crystal.py"] + [tool.ruff] line-length = 100 target-version = "py311" +[tool.ruff.lint] +# Correctness rules, not house style: pyflakes (unused names, duplicate dict +# keys, undefined names), syntax errors, comparison-with-itself, and loop +# variables captured by closures. +select = ["F", "E9", "PLR0124", "B023"] + [tool.mypy] python_version = "3.11" warn_unused_ignores = true disallow_untyped_defs = true ignore_missing_imports = true -[tool.pytest.ini_options] -testpaths = ["tests"] +# No pytest config: tests/ holds executable checks that take a live server and +# report by exit code, not pytest test functions. Configuring pytest over a +# directory with no tests gives a step that always passes while running +# nothing. diff --git a/queries/dimensionality.py b/queries/dimensionality.py deleted file mode 100644 index 7c88ef5..0000000 --- a/queries/dimensionality.py +++ /dev/null @@ -1,225 +0,0 @@ -"""Hydrogen-bond network dimensionality. - -What a crystallographer gets from this --------------------------------------- -Seed on one atom, walk the hydrogen-bond network outward, and count how many -atoms are reachable at each depth. The *shape of the growth curve* classifies -the packing: - - saturates after 2-4 atoms isolated dimer (0D) - roughly linear in depth chain (1D) - roughly quadratic sheet (2D) - roughly cubic framework (3D) - -That classification is not a curiosity. It predicts crystal morphology (chains -give needles, sheets give plates), cleavage and tabletting behaviour, and it -correlates with solubility and humidity stability. Deciding it today means -recomputing the packing geometrically for every candidate structure and then -running a connectivity analysis over the result. Here it is a traversal. - -Why this is the graph argument ------------------------------- -The hydrogen bonds themselves are the point. COD and the CSD store atoms and -coordinates; the contact network is *derived* every time anyone asks. We -computed it once at ingest, with the generating symmetry operation on every -edge, so the packing question becomes reachability rather than geometry. - -Dialect notes -------------- -TuringDB rejects an edge-type filter on a variable-length path -(`Edge type filters are not supported with variable-length paths`), so this runs -against a contact-only projection of the graph in which an untyped quantifier -*is* a hydrogen-bond traversal. The quantifier is TuringDB's postfix form -`-[e]->{k,k}`, not Neo4j's `-[e*k..k]->`. - -There is also no `length(path)` and no grouping aggregate, so the growth curve is -built by issuing one exact-depth query per depth and collecting the counts -client-side. Each query reports its own server-side execution time. -""" - -from __future__ import annotations - -import argparse -import json -import statistics -import sys -from dataclasses import dataclass, field -from typing import Sequence - -from turingdb import TuringDB - -DEFAULT_HOST = "http://localhost:6691" -DEFAULT_GRAPH = "cod_contacts" -MAX_DEPTH = 8 - - -@dataclass -class Curve: - """Reachability growth from one seed atom.""" - - seed_uid: str - cod_id: int | None - reached: list[int] = field(default_factory=list) # cumulative, per depth - ms: list[float] = field(default_factory=list) - - @property - def total(self) -> int: - return self.reached[-1] if self.reached else 0 - - def classify(self) -> tuple[str, str]: - """Classify the packing motif from the growth curve. - - Deliberately conservative and explained rather than fitted: we compare - successive increments, because that is what distinguishes a chain (an - increment that stays roughly constant) from a sheet (one that grows - roughly linearly) from a framework (one that keeps accelerating). - """ - r = self.reached - if not r or r[-1] < 2: - return "isolated", "no hydrogen-bonded neighbours found" - # per-depth increments - inc = [r[0]] + [r[i] - r[i - 1] for i in range(1, len(r))] - tail = [x for x in inc[2:] if x > 0] - if r[-1] <= 4 or not tail: - return "0D dimer / finite cluster", ( - f"growth stops at {r[-1]} atoms - an isolated motif, " - "typically more soluble" - ) - # is the increment still growing late in the walk? - early = statistics.mean(inc[2:4]) if len(inc) >= 4 else inc[-1] - late = statistics.mean(inc[-2:]) - ratio = late / early if early else 1.0 - if ratio < 1.35: - return "1D chain", ( - "reachability grows by a roughly constant amount per hop, " - "which is a chain - often needle morphology" - ) - if ratio < 2.5: - return "2D sheet", ( - "the per-hop increment itself grows roughly linearly, which is " - "a layer - plate morphology, easy cleavage" - ) - return "3D framework", ( - "the increment keeps accelerating, which is a three-dimensional " - "network - typically harder and less soluble" - ) - - -CYPHER_AT_DEPTH = ( - # Untyped quantifier: in cod_contacts the only edges are hydrogen bonds, - # so this IS a hydrogen-bond traversal. TuringDB does not accept - # -[e:CONTACT]->{{k,k}}. - "MATCH (a:Atom {{uid:'{uid}'}})-[e]->{{{k},{k}}}(b:Atom) RETURN count(b)" -) - - -def growth_curve( - client: TuringDB, seed_uid: str, max_depth: int = MAX_DEPTH -) -> Curve: - """Cumulative reachable-atom count at each depth from one seed.""" - curve = Curve(seed_uid=seed_uid, cod_id=None) - seen_total = 0 - for k in range(1, max_depth + 1): - q = CYPHER_AT_DEPTH.format(uid=seed_uid, k=k) - try: - n = int(client.query(q).iloc[0, 0]) - ms = client.get_query_exec_time() or 0.0 - except Exception: # noqa: BLE001 - break - # exact-depth counts are walks, not distinct atoms; accumulate the - # monotone envelope so the curve is interpretable as "reached by depth k" - seen_total = max(seen_total, n) - curve.reached.append(seen_total) - curve.ms.append(round(ms, 3)) - return curve - - -def pick_seeds(client: TuringDB, n: int) -> list[tuple[str, int]]: - """Atoms that actually participate in a hydrogen bond.""" - df = client.query( - "MATCH (a:Atom)-[e]->(b:Atom) RETURN a.uid, a.cod_id LIMIT " + str(n * 40) - ) - out: list[tuple[str, int]] = [] - seen: set[str] = set() - per_structure: set[int] = set() - for uid, cod in zip(df.iloc[:, 0], df.iloc[:, 1]): - uid = str(uid) - try: - cod_i = int(cod) - except Exception: # noqa: BLE001 - continue - # one seed per structure keeps the sample from being 40 atoms of the - # same crystal - if uid in seen or cod_i in per_structure: - continue - seen.add(uid) - per_structure.add(cod_i) - out.append((uid, cod_i)) - if len(out) >= n: - break - return out - - -def main(argv: Sequence[str] | None = None) -> int: - ap = argparse.ArgumentParser(description=__doc__) - ap.add_argument("--host", default=DEFAULT_HOST) - ap.add_argument("--graph", default=DEFAULT_GRAPH) - ap.add_argument("--seeds", type=int, default=12) - ap.add_argument("--depth", type=int, default=MAX_DEPTH) - ap.add_argument("--uid", default=None, help="run a single named seed atom") - ap.add_argument("--json", action="store_true") - args = ap.parse_args(argv) - - client = TuringDB(host=args.host) - try: - client.query(f"LOAD GRAPH {args.graph}") - except Exception: # noqa: BLE001 - pass - client.set_graph(args.graph) - - n_atoms = int(client.query("MATCH (n:Atom) RETURN count(n)").iloc[0, 0]) - n_edges = int(client.query("MATCH ()-[r]->() RETURN count(r)").iloc[0, 0]) - print(f"graph {args.graph}: {n_atoms:,} atoms, {n_edges:,} hydrogen bonds\n") - - seeds = ([(args.uid, None)] if args.uid - else pick_seeds(client, args.seeds)) - if not seeds: - print("no hydrogen-bonded atoms found") - return 1 - - curves: list[Curve] = [] - header = " " + "".join(f"{d:>8}" for d in range(1, args.depth + 1)) - print(f"{'seed atom':<26}{'COD':>9} {header} motif") - print("-" * (26 + 9 + 8 * args.depth + 24)) - for uid, cod in seeds: - c = growth_curve(client, uid, args.depth) - c.cod_id = cod - curves.append(c) - motif, _why = c.classify() - cells = "".join(f"{v:>8,}" for v in c.reached) - pad = " " * (8 * (args.depth - len(c.reached))) - print(f"{uid:<26}{str(cod or '-'):>9} {cells}{pad} {motif}") - - all_ms = [m for c in curves for m in c.ms] - if all_ms: - all_ms.sort() - print(f"\nserver-side latency over {len(all_ms)} traversals: " - f"p50 {all_ms[len(all_ms) // 2]:.2f} ms, " - f"p95 {all_ms[int(len(all_ms) * 0.95)]:.2f} ms, " - f"max {all_ms[-1]:.2f} ms") - - print("\ninterpretation:") - for c in curves[:6]: - motif, why = c.classify() - print(f" {c.seed_uid:<26} {motif:<26} {why}") - - if args.json: - print(json.dumps([{ - "seed": c.seed_uid, "cod_id": c.cod_id, "reached": c.reached, - "ms": c.ms, "motif": c.classify()[0], - } for c in curves], indent=2)) - return 0 - - -if __name__ == "__main__": - sys.exit(main()) diff --git a/queries/dimensionality_net.py b/queries/dimensionality_net.py index 752dce2..1b679b8 100644 --- a/queries/dimensionality_net.py +++ b/queries/dimensionality_net.py @@ -44,139 +44,12 @@ from pathlib import Path from typing import Sequence -import gemmi -import numpy as np - -LABELS = {0: "0D finite motif", 1: "1D chain", 2: "2D sheet", 3: "3D framework"} -CONSEQUENCE = { - 0: "isolated dimers/clusters - typically the most soluble", - 1: "chains - often needle morphology, anisotropic mechanical response", - 2: "layers - plate morphology, easy cleavage, good tabletting", - 3: "three-dimensional network - harder, less soluble, humidity stable", -} - - -def _compose(g: gemmi.Op, s: gemmi.Op) -> tuple[np.ndarray, np.ndarray]: - """Compose g after s, KEEPING the lattice part of the translation. - - gemmi's own `g * s` reduces the translation modulo one unit cell, which - silently discards exactly the information this calculation depends on: - - S = -x,-y+1,-z tran [0, 24, 0] - g[0] * S = -x,-y,-z tran [0, 0, 0] <- the +1 in y is gone - - With every lattice translation flattened to zero, every cycle closes with a - zero period and every structure looks 0D. So we do the composition in exact - integer arithmetic instead. gemmi scales both rot and tran by Op.DEN, so - with R = rot/DEN an integer rotation matrix: - - rot(g.s) = R_g @ rot_s (still DEN-scaled) - tran(g.s) = R_g @ tran_s + tran_g (still DEN-scaled, NOT reduced) - """ - den = gemmi.Op.DEN - rg = np.array(g.rot, dtype=np.int64) // den - rot = rg @ np.array(s.rot, dtype=np.int64) - tran = rg @ np.array(s.tran, dtype=np.int64) + np.array(g.tran, dtype=np.int64) - return rot, tran - - -def _decompose( - rot: np.ndarray, tran: np.ndarray, ops: list[gemmi.Op] -) -> tuple[int, tuple[int, int, int]] | None: - """Write (rot, tran) as ops[k] followed by an integer lattice translation.""" - den = gemmi.Op.DEN - for k, op in enumerate(ops): - if not np.array_equal(np.array(op.rot, dtype=np.int64), rot): - continue - d = tran - np.array(op.tran, dtype=np.int64) - if np.all(d % den == 0): - q = d // den - return k, (int(q[0]), int(q[1]), int(q[2])) - return None - - -def _op_from_code(code: str, ops: list[gemmi.Op]) -> gemmi.Op: - idx_s, tr_s = code.split("_") - base = ops[int(idx_s) - 1] - op = gemmi.Op(base.triplet()) - den = gemmi.Op.DEN - sh = [int(c) - 5 for c in tr_s] - op.tran = (op.tran[0] + sh[0] * den, - op.tran[1] + sh[1] * den, - op.tran[2] + sh[2] * den) - return op - - -def net_dimensionality( - spacegroup_number: int, - contacts: list[tuple[str, str, str]], - uid_to_site: dict[str, int], -) -> tuple[int, int, int]: - """Return (dimensionality, n_nodes_in_largest_component, n_components). - - `contacts` is [(a_uid, b_uid, symop_code)]; a_uid sits at the identity and - b_uid is generated by symop_code. - """ - sg = gemmi.find_spacegroup_by_number(spacegroup_number) - if sg is None or not contacts: - return 0, 0, 0 - ops = list(sg.operations()) - n_ops = len(ops) - - # Unit-cell nodes are (asymmetric site, symmetry operation). Applying every - # operation g to a stored edge gives the whole orbit of that contact. - adj: dict[tuple[int, int], list[tuple[tuple[int, int], np.ndarray]]] = defaultdict(list) - for a_uid, b_uid, code in contacts: - a = uid_to_site.get(a_uid) - b = uid_to_site.get(b_uid) - if a is None or b is None: - continue - try: - s_op = _op_from_code(code, ops) - except Exception: # noqa: BLE001 - continue - for g_idx, g in enumerate(ops): - c_rot, c_tran = _compose(g, s_op) - got = _decompose(c_rot, c_tran, ops) - if got is None: - continue - k, tau = got - v = np.array(tau, dtype=np.int64) - adj[(a, g_idx)].append(((b, k), v)) - adj[(b, k)].append(((a, g_idx), -v)) - - if not adj: - return 0, 0, 0 - - # spanning forest with accumulated lattice potentials; every non-tree edge - # closes a cycle whose net translation is a period of the net - seen: dict[tuple[int, int], np.ndarray] = {} - cycles: list[np.ndarray] = [] - comp_sizes: list[int] = [] - for start in list(adj): - if start in seen: - continue - seen[start] = np.zeros(3, dtype=np.int64) - stack = [start] - size = 0 - while stack: - u = stack.pop() - size += 1 - for v, t in adj.get(u, ()): - if v not in seen: - seen[v] = seen[u] + t - stack.append(v) - else: - cyc = seen[u] + t - seen[v] - if cyc.any(): - cycles.append(cyc) - comp_sizes.append(size) - - if not cycles: - return 0, max(comp_sizes, default=0), len(comp_sizes) - rank = int(np.linalg.matrix_rank(np.array(cycles, dtype=np.float64))) - return min(rank, 3), max(comp_sizes, default=0), len(comp_sizes) - +# The calculation itself lives in ingest/netdim.py so that build_graph.py can +# stamp net_dim onto every Structure node at ingest time and the UI can query it +# directly. This module is the CLI over it: corpus-wide distribution, plus the +# strong-vs-weak comparison that justifies admitting C-H...A at all. +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from ingest.netdim import CONSEQUENCE, LABELS, net_dimensionality # noqa: E402 def main(argv: Sequence[str] | None = None) -> int: ap = argparse.ArgumentParser(description=__doc__) @@ -185,6 +58,12 @@ def main(argv: Sequence[str] | None = None) -> int: ap.add_argument("--limit", type=int, default=None) ap.add_argument("--include-inferred", action="store_true", help="include hydrogen bonds with no located hydrogen") + ap.add_argument("--kinds", default="hbond", + help="comma-separated contact kinds to traverse. " + "'hbond' is the strong population; add 'hbond_weak' " + "to admit C-H...A, which is the comparison that " + "decides whether the 0D share is real or an artefact " + "of the criteria") ap.add_argument("--out", type=Path, default=Path("data/manifest/dimensionality.json")) args = ap.parse_args(argv) @@ -193,6 +72,8 @@ def main(argv: Sequence[str] | None = None) -> int: print(f"missing {args.cache}", file=sys.stderr) return 1 + kinds = {k.strip() for k in args.kinds.split(",") if k.strip()} + dist: Counter[int] = Counter() examples: dict[int, list[tuple[int, str, int]]] = defaultdict(list) rows: list[dict] = [] @@ -212,7 +93,7 @@ def main(argv: Sequence[str] | None = None) -> int: contacts = [ (c["a"], c["b"], c["symop"]) for c in res["contacts"] - if c["kind"] == "hbond" + if c["kind"] in kinds and (args.include_inferred or not c["h_inferred"]) ] if not contacts: @@ -231,7 +112,8 @@ def main(argv: Sequence[str] | None = None) -> int: (res["cod_id"], res["spacegroup_hm"], len(contacts)) ) - print(f"structures with at least one located-H hydrogen bond: {n:,}") + print(f"contact kinds traversed: {', '.join(sorted(kinds))}") + print(f"structures with at least one qualifying hydrogen bond: {n:,}") print(f"(inferred-H bonds {'INCLUDED' if args.include_inferred else 'excluded'} " f"- a motif claim resting on unrefined hydrogens would not survive review)\n") print(f"{'dimensionality':<20}{'structures':>12}{'share':>9} consequence") @@ -249,6 +131,7 @@ def main(argv: Sequence[str] | None = None) -> int: args.out.write_text(json.dumps( {"counts": {LABELS[k]: v for k, v in sorted(dist.items())}, "n_structures": n, "include_inferred": args.include_inferred, + "kinds": sorted(kinds), "per_structure": rows}, indent=2) + "\n") print(f"\nwrote {args.out}") return 0 diff --git a/run.sh b/run.sh index 09b4456..d21eb13 100755 --- a/run.sh +++ b/run.sh @@ -1,11 +1,14 @@ #!/usr/bin/env bash -# One-command setup for the CCDC crystallography demo. +# One-command setup. # -# ./run.sh start TuringDB + the visualizer (assumes data is built) +# ./run.sh start TuringDB + the browser studio on existing data # ./run.sh --ingest also acquire COD and build the graphs from scratch # -# The ingest path takes a while and needs network access to +# The ingest path takes roughly 25 minutes and needs network access to # crystallography.net (rsync for CIFs, MySQL for metadata). +# +# Once running: http://localhost:8087 for the studio, `crystal facts` for the +# command line. See TUTORIAL.md. set -euo pipefail cd "$(dirname "$0")" @@ -38,21 +41,65 @@ if [[ "${1:-}" == "--ingest" ]]; then echo "==> acquiring COD (CIFs by rsync, metadata from the public MySQL mirror)" uv run python -m ingest.download --dataset slice echo "==> building the graph" - uv run python -m ingest.build_graph --dataset slice --limit 10000 - uv run python -m ingest.project_contacts + # --out is explicit: build_graph would otherwise emit data/jsonl/slice.jsonl + # while load_turingdb below reads cod_slice.jsonl, and the mismatch only + # shows up as a missing file at the end of a long ingest. + uv run python -m ingest.build_graph --dataset slice \ + --out data/jsonl/cod_slice.jsonl + uv run python -m ingest.project_contacts \ + --in data/jsonl/cod_slice.jsonl --out data/jsonl/cod_contacts_v2.jsonl echo "==> loading into TuringDB" TURING_HOST="http://localhost:$TURING_PORT" \ uv run python -m ingest.load_turingdb \ --jsonl data/jsonl/cod_slice.jsonl --graph cod_slice_v2 \ --host "http://localhost:$TURING_PORT" --turing-dir "$TURING_DIR" + TURING_HOST="http://localhost:$TURING_PORT" \ + uv run python -m ingest.load_turingdb \ + --jsonl data/jsonl/cod_contacts_v2.jsonl --graph cod_contacts_v2 \ + --host "http://localhost:$TURING_PORT" --turing-dir "$TURING_DIR" echo "==> building the chronologically-committed corpus (one commit per year)" uv run python -m ingest.build_versioned \ --host "http://localhost:$TURING_PORT" --per-year 0 --fresh fi +# A graph that exists on disk is NOT resident after a server start, and an +# unloaded graph answers GRAPH_NOT_FOUND rather than returning nothing -- which +# reads as "the demo is broken" rather than "load it first". Load whatever is +# there, and say plainly when there is nothing. +echo "==> loading graphs" +uv run python - "$TURING_PORT" <<'PY' +import sys +from turingdb import TuringDB + +host = f"http://localhost:{sys.argv[1]}" +c = TuringDB(host=host) +try: + available = set(c.list_available_graphs()["graphName"]) +except Exception as exc: + print(f" cannot reach TuringDB at {host}: {exc}") + raise SystemExit(1) + +wanted = [g for g in ("cod_slice_v2", "cod_contacts_v2", + "cod_versioned", "cod_versions") if g in available] +if not wanted: + print(" no crystallography graphs on this server yet.") + print(" Run ./run.sh --ingest to acquire COD and build them") + print(" (~2.8 GB of CIFs, roughly 25 minutes).") + raise SystemExit(0) + +for g in wanted: + try: + c.load_graph(g) + print(f" loaded {g}") + except Exception as exc: + # already loaded is the common, harmless case + print(f" {g}: {str(exc)[:80]}") +PY + echo "==> frontend" cd visualizer [ -d node_modules ] || npm install npm run build -echo "==> serving on :$UI_PORT (proxying /api -> :$TURING_PORT)" +echo "==> studio on http://localhost:$UI_PORT (proxying /api -> :$TURING_PORT)" +echo "==> command line: TURING_HOST=http://localhost:$TURING_PORT uv run python crystal.py facts" TURING_FRONTEND_PORT="$UI_PORT" TURING_API_PORT="$TURING_PORT" node server.js diff --git a/scripts/tunnel.sh b/scripts/tunnel.sh deleted file mode 100755 index 8f2da3f..0000000 --- a/scripts/tunnel.sh +++ /dev/null @@ -1,13 +0,0 @@ -#!/bin/bash -# Self-healing tunnel for the CCDC crystallography demo. -# 8087 -> visualizer (proxies /api to TuringDB on 6691) -# 6691 -> TuringDB REST, for running Cypher directly from the laptop -# macOS has no setsid, so launch with nohup. -while true; do - ssh -N -o ExitOnForwardFailure=yes -o ServerAliveInterval=30 \ - -L 8087:localhost:8087 \ - -L 6691:localhost:6691 \ - ubuntu@hetz-adam - echo "$(date -Is) tunnel dropped, reconnecting in 3s" - sleep 3 -done diff --git a/tests/check_cli_queries.py b/tests/check_cli_queries.py new file mode 100644 index 0000000..060cf7d --- /dev/null +++ b/tests/check_cli_queries.py @@ -0,0 +1,52 @@ +"""Execute every named query the CLI ships, against a live server. + +Same purpose as tests/check_frontend_queries.py, for the other surface. A +query that no longer parses -- because a property was renamed, or a graph was +rebuilt without it -- should fail here rather than in front of someone. + + python tests/check_cli_queries.py + python tests/check_cli_queries.py --host http://localhost:6691 +""" + +from __future__ import annotations + +import argparse +import sys +import time +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) + +import crystal # noqa: E402 + + +def main(argv: list[str] | None = None) -> int: + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--host", default=crystal.DEFAULT_HOST) + ap.add_argument("--graph", default=crystal.DEFAULT_GRAPH) + a = ap.parse_args(argv) + + failed: list[tuple[str, str]] = [] + print(f"{'query':<18}{'rows':>10}{'ms':>9}") + print("-" * 37) + for name, (_desc, cypher, _tally) in crystal.LIBRARY.items(): + t0 = time.perf_counter() + try: + _, rows = crystal.query(cypher, graph=a.graph, host=a.host) + except Exception as exc: # noqa: BLE001 + failed.append((name, str(exc)[:160])) + print(f"{name:<18}{'ERROR':>10}") + continue + print(f"{name:<18}{len(rows):>10,}{(time.perf_counter() - t0) * 1000:>9.0f}") + + if failed: + print(f"\n{len(failed)} FAILING:") + for name, err in failed: + print(f" {name}: {err}") + return 1 + print(f"\nall {len(crystal.LIBRARY)} named queries execute cleanly") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/check_frontend_queries.py b/tests/check_frontend_queries.py new file mode 100644 index 0000000..f406d96 --- /dev/null +++ b/tests/check_frontend_queries.py @@ -0,0 +1,328 @@ +"""Execute every Cypher string the frontend ships, against a live server. + +The failure mode this exists to catch is specific to this dialect and this app: +a projection of an absent property is an `ANALYZE_ERROR`, not an empty column, +and the UI wraps several queries in a `.catch(() => [])`. So a query that broke +because a property was renamed, or because a graph was rebuilt without it, +renders as **an empty canvas rather than an error** -- and an empty canvas reads +as "no data" to everyone in the room. + +That has already bitten this family of demos more than once: `type(r)` became a +PARSE_ERROR in 1.36 and a swallowed error drew nodes with no links; here, +`IN_COMPONENT` and `net_dim` exist only in the v2 corpus, so four of the shipped +investigations would silently return nothing against an older graph. + +The check is deliberately dumb: pull every template literal and string that +looks like Cypher out of the crystal module, run it, and require a non-error +response. Negative-control it by injecting a known-bad query before trusting a +green run. + +Usage: + .venv/bin/python tests/check_frontend_queries.py --graph cod_slice_v2 +""" + +from __future__ import annotations + +import argparse +import json +import re +import sys +import urllib.request +from pathlib import Path + +CRYSTAL = Path("visualizer/src/components/viewer/crystal") + +#: A string is treated as Cypher only if it STARTS a query. `RETURN` and +#: `WHERE` are deliberately absent: they occur mid-string in concatenated +#: fragments and in prose, and including them made this check report its own +#: parsing as failures rather than the app's. +STARTS = ("MATCH ", "MATCH(", "CALL ", "LOAD ") + +#: Anything containing one of these is JavaScript that merely begins with a +#: Cypher keyword -- a concatenation, a JSX block, an arrow function. +#: NOTE: "${" is deliberately NOT here. Interpolations are resolved by SUBS +#: below; rejecting them at this stage discards every non-trivial query and +#: leaves the check silently covering only the literal ones. +NOT_CYPHER = ('")', "\")", "=>", "' + +#: Placeholders the UI substitutes at runtime. Replaced with something valid so +#: the query shape can still be executed. +SUBS = { + r"\$\{PACKING_SEED_COD\}": "2229029", + r"\$\{cod\}": "2229029", + r"\$\{hm\}": "P -1", + r"\$\{t\}": "carboxylic_acid", + r"\$\{id\}": "0", + r"\$\{proj\}": "n.cod_id", + r"\$\{year\}": "2011", + r"\$\{k\}": "2", + # DIMER_COUNT and base are query PREFIXES concatenated with a suffix at the + # call site; substituting their text is what makes those queries runnable. + r"\$\{DIMER_COUNT\}": ( + "MATCH (f1:Fragment)<-[:IN_FRAGMENT]-(a1:Atom)" + "-[h:CONTACT]->(a2:Atom)-[:IN_FRAGMENT]->(f2:Fragment) " + "WHERE f1.fragment_type = 'carboxylic_acid' " + "AND f2.fragment_type = 'carboxylic_acid' " + "AND h.kind = 'hbond' AND h.h_inferred = false AND h.is_involution = "), + r"\$\{base\}": ( + "MATCH (g1:SpaceGroup)<-[:IN_SPACE_GROUP]-(s1:Structure)" + "-[r1:CONTAINS_COMPONENT]->(c:Component)" + "<-[r2:CONTAINS_COMPONENT]-(s2:Structure)-[:IN_SPACE_GROUP]->(g2:SpaceGroup) " + "WHERE r1.role = 'principal' AND r2.role = 'principal' " + "AND c.is_solvent = false AND c.has_inchikey = true "), +} + + +#: The investigation queries: `cypher: `...`` fields. These are the queries +#: that paint the canvas and the ones a schema change breaks first, and the +#: field is well delimited, so they can be recovered exactly rather than +#: guessed at. +_CYPHER_FIELD = re.compile(r"cypher:\s*`([^`]*)`", re.S) + +#: Single-literal queries passed directly to a helper, e.g. +#: `scalar(g, "MATCH ... RETURN count(n)")`. Multi-literal concatenations are +#: NOT recovered -- see the coverage assertion below, which exists precisely so +#: that gap is reported rather than hidden. +_CALL_QUERY = re.compile( + r"(?:scalar|scalarAt|runCypher|queryAt|loadSubgraph|q)\(\s*" + r"(?:[A-Za-z_][\w.]*\s*,\s*)?" + r"(?:\"((?:[^\"\\]|\\.)*)\"|'((?:[^'\\]|\\.)*)'|`([^`]*)`)\s*\)", + re.S) + +#: Tokens that MUST appear somewhere in the extracted set. Each is a schema +#: element the studio depends on and that has no other coverage: if a rebuild +#: drops one, the panels render empty rather than failing, which is the exact +#: class of bug this file exists to catch. An empty extraction therefore fails +#: loudly instead of reporting success over nothing. +#: Schema elements the STUDIO depends on. IN_COMPONENT is deliberately absent: +#: the studio's motif queries reach functional groups through IN_FRAGMENT, which +#: is atom-level and therefore correct, so nothing here traverses atom -> +#: molecule. That edge is covered by tests/check_cli_queries.py instead. +REQUIRED_TOKENS = ("CONTACT", "HAS_FRAGMENT", "IN_FRAGMENT", + "net_dim", "is_involution", "Snapshot", "IN_SPACE_GROUP") + + +def extract(paths: list[Path]) -> list[tuple[str, str]]: + out: list[tuple[str, str]] = [] + for p in paths: + text = p.read_text() + cands: list[tuple[int, str]] = [] + for m in _CYPHER_FIELD.finditer(text): + cands.append((m.start(), m.group(1))) + for m in _CALL_QUERY.finditer(text): + cands.append((m.start(), m.group(1) or m.group(2) or m.group(3) or "")) + for pos, raw in cands: + body = raw.replace("\\'", "'").replace('\\"', '"').strip() + if not body.upper().startswith(STARTS): + continue + if any(tok in body for tok in NOT_CYPHER): + continue + for pat, rep in SUBS.items(): + body = re.sub(pat, rep, body) + if "${" in body: + continue + line = text[:pos].count("\n") + 1 + out.append((f"{p.name}:{line}", " ".join(body.split()))) + return out + + +def _count_entries(text: str, i: int) -> int: + """Count top-level entries of an array literal whose `[` is just before `i`. + + Written as an explicit scanner rather than a regex because the entries are + `scalar(g, `...${x}...`)` calls: they nest parentheses, contain commas + inside both ordinary strings and template literals, and template literals + contain `${ }` holes that themselves nest. Anything less careful + miscounts, which would make this check worse than useless -- it would + report mismatches that are not real and train the reader to ignore it. + """ + depth = 1 + entries = 0 + saw_content = False + last_meaningful = "" + n = len(text) + while i < n and depth > 0: + ch = text[i] + if ch in "'\"": # ordinary string + q, i = ch, i + 1 + while i < n and text[i] != q: + i += 2 if text[i] == "\\" else 1 + saw_content, last_meaningful = True, "x" + elif ch == "`": # template literal, may nest ${} + i += 1 + hole = 0 + while i < n: + c = text[i] + if c == "\\": + i += 2 + continue + if c == "$" and i + 1 < n and text[i + 1] == "{": + hole += 1 + i += 2 + continue + if c == "}" and hole: + hole -= 1 + elif c == "`" and not hole: + break + i += 1 + saw_content, last_meaningful = True, "x" + elif ch == "/" and i + 1 < n and text[i + 1] == "/": + while i < n and text[i] != "\n": # line comment + i += 1 + continue + elif ch in "([{": + depth += 1 + saw_content, last_meaningful = True, "x" + elif ch in ")]}": + depth -= 1 + if depth == 0: + break + saw_content, last_meaningful = True, "x" + elif ch == "," and depth == 1: + entries += 1 + last_meaningful = "," + elif not ch.isspace(): + saw_content, last_meaningful = True, "x" + i += 1 + if not saw_content: + return 0 + # a trailing comma before `]` means entries == item count, else count - 1 + return entries if last_meaningful == "," else entries + 1 + + +def check_destructure_arity(paths: list[Path]) -> list[str]: + """`const [a, b, c] = await Promise.all([...])` must name exactly as many + variables as the array has entries. + + This is not a style check. The panels read their numbers POSITIONALLY, so + inserting one query into the wrong Promise.all shifts every variable after + it: a count of structures lands in the variable that prints as + "carboxylate components", and the panel renders a plausible, wrong number + with no error at all. The only symptom may be a crash in a *different* + panel whose array has gone one short. + """ + bad: list[str] = [] + pat = re.compile(r"const\s*\[([^\]]*)\]\s*=\s*await\s+Promise\.all\(\[", re.S) + for p in paths: + text = p.read_text() + for m in pat.finditer(text): + names = [v.strip() for v in m.group(1).split(",") if v.strip()] + count = _count_entries(text, m.end()) + if count != len(names): + line = text[: m.start()].count("\n") + 1 + bad.append( + f"{p.name}:{line} destructures {len(names)} " + f"({', '.join(names)}) from a Promise.all of {count}") + return bad + + +def run(host: str, graph: str, cypher: str) -> str: + req = urllib.request.Request( + f"{host}/query?graph={graph}", + data=cypher.encode(), + headers={"Content-Type": "application/json"}, + method="POST", + ) + try: + with urllib.request.urlopen(req, timeout=300) as r: + payload = json.loads(r.read()) + except Exception as exc: # noqa: BLE001 + return f"{type(exc).__name__}: {exc}"[:120] + # Failures arrive as HTTP 200 with a non-null top-level `error`. + if payload.get("error"): + return f"{payload['error']}: {payload.get('error_details', '')}"[:140] + return "" + + +def main(argv: list[str] | None = None) -> int: + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--host", default="http://localhost:6691") + ap.add_argument("--graph", default="cod_slice_v2") + ap.add_argument("--ledger-graph", default="cod_versions") + ap.add_argument("--versioned-graph", default="cod_versioned") + ap.add_argument("--negative-control", action="store_true", + help="inject a known-bad query and require that it FAILS") + a = ap.parse_args(argv) + + files = sorted(CRYSTAL.glob("*.ts")) + sorted(CRYSTAL.glob("*.tsx")) + if not files: + print(f"no frontend sources under {CRYSTAL}", file=sys.stderr) + return 1 + arity = check_destructure_arity(files) + if arity: + print("DESTRUCTURE / Promise.all ARITY MISMATCH -- panels read these " + "positionally, so this renders wrong numbers, not an error:") + for a_ in arity: + print(f" {a_}") + return 1 + + queries = extract(files) + + if a.negative_control: + queries.append(("NEGATIVE-CONTROL", + "MATCH (n:Structure) RETURN n.property_that_does_not_exist")) + + # Snapshot queries belong to the ledger graph, Snapshot-free ones to the + # corpus. Route by the labels they mention rather than by guessing. + def graph_for(cypher: str) -> str: + if "Snapshot" in cypher: + return a.ledger_graph + return a.graph + + bad: list[tuple[str, str, str]] = [] + missing: dict[str, int] = {} + seen: set[str] = set() + n = 0 + for where, cypher in queries: + if cypher in seen: + continue + seen.add(cypher) + n += 1 + target = graph_for(cypher) + err = run(a.host, target, cypher) + # A graph that is absent entirely (mid-rebuild, or never built on this + # server) is reported separately: it is a deployment fact, not a broken + # query, and conflating the two makes the check cry wolf. + if err.startswith("GRAPH_NOT_FOUND"): + missing.setdefault(target, 0) + missing[target] += 1 + continue + if err: + bad.append((where, cypher, err)) + + print(f"{n} distinct Cypher strings from {len(files)} frontend files") + uncovered = [t for t in REQUIRED_TOKENS + if not any(t in c for _, c in queries)] + if uncovered: + print("\nNOT COVERED by any extracted query: " + ", ".join(uncovered)) + print("Those schema elements are used by the studio but no query " + "reaching them was recovered, so this run proves nothing about " + "them. Fix the extractor rather than trusting the result.") + if a.negative_control: + ctrl = [b for b in bad if b[0] == "NEGATIVE-CONTROL"] + if not ctrl: + print("NEGATIVE CONTROL PASSED WHEN IT SHOULD HAVE FAILED -- " + "this check is not actually executing anything", file=sys.stderr) + return 2 + print("negative control failed as required (the check works)") + bad = [b for b in bad if b[0] != "NEGATIVE-CONTROL"] + + for g, k in sorted(missing.items()): + print(f" NOTE: {k} queries skipped -- graph '{g}' is not loaded") + if not bad and not uncovered: + print(f"all {n} extracted queries execute cleanly, and every required " + f"schema element is covered") + return 0 + if not bad: + print(f"\nthe {n} extracted queries execute cleanly, but coverage is " + f"incomplete (above)") + return 1 + print(f"\n{len(bad)} FAILING:") + for where, cypher, err in bad: + print(f"\n {where}\n {cypher[:170]}\n -> {err}") + return 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/check_slow_structure.py b/tests/check_slow_structure.py deleted file mode 100644 index 83aa04f..0000000 --- a/tests/check_slow_structure.py +++ /dev/null @@ -1,55 +0,0 @@ -"""Regression check for the structures that used to stall the ingest. - -COD 4504526 (147 sites, P2(1)/c) has an 85-atom organic component whose -bond-order perception never returned with RDKit's default unlimited iteration -budget. Kept as a named test so the fix does not silently regress. -""" - -from __future__ import annotations - -import sys -import time -from pathlib import Path - -from ingest.build_graph import process_structure -from ingest.download import cod_cif_relpath - -#: (cod_id, generous wall-clock budget in seconds) -CASES: list[tuple[int, float]] = [ - (4504526, 30.0), -] - - -def main() -> int: - ok = True - for cod_id, budget in CASES: - path = Path("data/cif") / cod_cif_relpath(cod_id) - if not path.exists(): - print(f" SKIP COD {cod_id} (CIF not on disk)") - continue - t0 = time.perf_counter() - res = process_structure((str(path), cod_id)) - dt = time.perf_counter() - t0 - if res is None: - print(f" FAIL COD {cod_id}: returned None") - ok = False - continue - skipped = res.get("skipped") - n_comp = len(res.get("components", [])) - n_con = len(res.get("contacts", [])) - n_ik = sum(1 for c in res.get("components", []) if c.get("inchikey")) - verdict = "PASS" if dt < budget else "FAIL" - if dt >= budget: - ok = False - print(f" {verdict} COD {cod_id}: {dt:.2f}s (budget {budget:.0f}s) " - f"components={n_comp} with_inchikey={n_ik} contacts={n_con} " - f"skipped={skipped}") - if res.get("stats", {}).get("timed_out"): - print(" note: contact search hit its time budget " - "(contacts_truncated will be set)") - print("\nRESULT:", "PASS" if ok else "FAIL") - return 0 if ok else 1 - - -if __name__ == "__main__": - sys.exit(main()) diff --git a/tests/corpus_facts.py b/tests/corpus_facts.py new file mode 100644 index 0000000..2d6f4d0 --- /dev/null +++ b/tests/corpus_facts.py @@ -0,0 +1,150 @@ +"""Re-derive every number the README and the UI quote, from the live graph. + +Written because a review found that three of the five headline figures in the +README were computed over the wrong population, and nothing in the repo would +have caught it: + +* "mean H...A 2.09 A at 160.5 deg" filtered on `has_angle`, which **includes the + halogen bonds** -- so it averaged H...A together with X...A. The located-H + hydrogen-bond value is 2.053 A / 160.30 deg. +* "centrosymmetric groups account for 74.0%" was the sum of five hand-picked + space groups. Asked of the data it is 82.4%. +* "33,930 intermolecular contacts" blends three populations that the demo's own + rules say must never be blended. + +The fix is not to correct the prose once -- it is to make the prose quote a +script, so the next rebuild cannot drift from it. Run this after every rebuild +and paste the block it prints. + +Usage: + .venv/bin/python tests/corpus_facts.py --graph cod_slice_v2 +""" + +from __future__ import annotations + +import argparse +import json +import sys +import urllib.request + + +def scalar(host: str, graph: str, cypher: str) -> float: + req = urllib.request.Request( + f"{host}/query?graph={graph}", + data=cypher.encode(), + headers={"Content-Type": "application/json"}, + method="POST", + ) + with urllib.request.urlopen(req, timeout=600) as r: + payload = json.loads(r.read()) + if payload.get("error"): + raise RuntimeError(f"{payload['error']}: {payload.get('error_details')}") + data = payload.get("data") or [[[None]]] + try: + v = data[0][0][0] + except (IndexError, TypeError): + return 0.0 + return float(v) if v is not None else 0.0 + + +def main(argv: list[str] | None = None) -> int: + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--host", default="http://localhost:6691") + ap.add_argument("--graph", default="cod_slice_v2") + ap.add_argument("--out", default="data/manifest/corpus_facts.json") + a = ap.parse_args(argv) + S = lambda c: scalar(a.host, a.graph, c) # noqa: E731 + + f: dict[str, float] = {} + f["structures"] = S("MATCH (s:Structure) RETURN count(s)") + f["atoms"] = S("MATCH (n:Atom) RETURN count(n)") + f["components"] = S("MATCH (n:Component) RETURN count(n)") + f["bonds"] = S("MATCH ()-[r:BONDED_TO]->() RETURN count(r)") + f["space_groups"] = S("MATCH (n:SpaceGroup) RETURN count(n)") + + # --- the three contact populations, never blended --------------------- + f["hbond_strong"] = S( + "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' " + "AND r.h_inferred = false RETURN count(r)") + f["hbond_inferred"] = S( + "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' " + "AND r.h_inferred = true RETURN count(r)") + f["hbond_weak"] = S( + "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond_weak' RETURN count(r)") + f["halogen"] = S( + "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'halogen' RETURN count(r)") + f["close_contact"] = S( + "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'close_contact' RETURN count(r)") + f["contacts_all"] = S("MATCH ()-[r:CONTACT]->() RETURN count(r)") + + # --- geometry, over the located-H HYDROGEN-BOND population only ------- + base = ("MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' " + "AND r.h_inferred = false RETURN ") + f["mean_HA"] = S(base + "avg(r.length)") + f["mean_angle"] = S(base + "avg(r.angle)") + wbase = "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond_weak' RETURN " + f["mean_HA_weak"] = S(wbase + "avg(r.length)") + f["mean_angle_weak"] = S(wbase + "avg(r.angle)") + + # --- symmetry, asked of the data rather than a hand list -------------- + f["centrosymmetric"] = S( + "MATCH (s:Structure)-[:IN_SPACE_GROUP]->(x:SpaceGroup) " + "WHERE x.is_centrosymmetric = true RETURN count(s)") + + # --- periodic-net dimensionality -------------------------------------- + for k in (0, 1, 2, 3): + f[f"net_dim_{k}"] = S( + f"MATCH (s:Structure) WHERE s.net_dim = {k} RETURN count(s)") + f[f"net_dim_weak_{k}"] = S( + f"MATCH (s:Structure) WHERE s.net_dim_weak = {k} RETURN count(s)") + f["net_scored"] = sum(f[f"net_dim_{k}"] for k in range(4)) + f["net_scored_weak"] = sum(f[f"net_dim_weak_{k}"] for k in range(4)) + + f["with_inchikey"] = S( + "MATCH (c:Component) WHERE c.has_inchikey = true RETURN count(c)") + f["acid_components"] = S( + "MATCH (c:Component)-[:HAS_FRAGMENT]->(f:Fragment) " + "WHERE f.fragment_type = 'carboxylic_acid' RETURN count(c)") + + pct = lambda x, y: 100.0 * x / y if y else 0.0 # noqa: E731 + print(f"\n=== {a.graph} ===\n") + print(f"Structures {f['structures']:>12,.0f}") + print(f"Atomic sites {f['atoms']:>12,.0f}") + print(f"Components {f['components']:>12,.0f} " + f"InChIKey {pct(f['with_inchikey'], f['components']):.1f}%") + print(f"Covalent bonds {f['bonds']:>12,.0f}") + print(f"Space groups {f['space_groups']:>12,.0f} " + f"centrosymmetric {pct(f['centrosymmetric'], f['structures']):.1f}% " + f"of structures") + print("\nContacts -- three populations, never blended:") + print(f" strong H-bond (located H) {f['hbond_strong']:>12,.0f} " + f"mean H...A {f['mean_HA']:.3f} A at {f['mean_angle']:.2f} deg") + print(f" inferred (no located H) {f['hbond_inferred']:>12,.0f} " + f"{pct(f['hbond_inferred'], f['hbond_strong'] + f['hbond_inferred']):.1f}% " + f"of the H-bond layer") + print(f" weak C-H...A {f['hbond_weak']:>12,.0f} " + f"mean H...A {f['mean_HA_weak']:.3f} A at {f['mean_angle_weak']:.2f} deg") + print(f" halogen bond {f['halogen']:>12,.0f}") + print(f" close_contact (excluded) {f['close_contact']:>12,.0f}") + print(f" ALL CONTACT edges {f['contacts_all']:>12,.0f}") + print("\nPeriodic-net dimensionality (h_inferred excluded):") + lab = {0: "0D finite motif", 1: "1D chain", 2: "2D sheet", 3: "3D framework"} + for k in (0, 1, 2, 3): + print(f" {lab[k]:<18} strong {f[f'net_dim_{k}']:>8,.0f} " + f"({pct(f[f'net_dim_{k}'], f['net_scored']):>5.1f}%)" + f" +weak {f[f'net_dim_weak_{k}']:>8,.0f} " + f"({pct(f[f'net_dim_weak_{k}'], f['net_scored_weak']):>5.1f}%)") + print(f" scored: {f['net_scored']:,.0f} strong / " + f"{f['net_scored_weak']:,.0f} with weak\n") + + from pathlib import Path + + out = Path(a.out) + out.parent.mkdir(parents=True, exist_ok=True) + out.write_text(json.dumps(f, indent=2) + "\n") + print(f"wrote {out}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/measure_perf.py b/tests/measure_perf.py new file mode 100644 index 0000000..e2818be --- /dev/null +++ b/tests/measure_perf.py @@ -0,0 +1,270 @@ +"""Performance numbers for the CCDC talk, measured rather than asserted. + +The Journal Club brief asks for "performance, scalability, ease of use" by name, +and the honest way to answer that is a table of absolute latencies at a stated +corpus size plus an extrapolation to CSD scale -- not a rigged head-to-head +against another engine. + +Reports four things: + +1. **Read latency** for every query the UI actually ships, plus the join-order + comparison, taken as the server's own reported execution time (median of N). +2. **Datapart count and bytes on disk** before and after `MERGE_DATAPARTS`. + A single-shot `LOAD JSONL` produces one datapart, so the merge has nothing + to do; the numbers are reported so that is visible rather than assumed. +3. **Property index effect** on the same queries. +4. **The CSD extrapolation**: bytes and elements per structure, scaled. + +Usage: + .venv/bin/python tests/measure_perf.py --graph cod_slice_v2 \\ + --turing-dir turing-data-137 --host http://localhost:6691 +""" + +from __future__ import annotations + +import argparse +import json +import statistics +import sys +import time +from pathlib import Path + +import urllib.request + +#: The CSD is roughly 1.3 M structures; this is what the database team will ask. +CSD_STRUCTURES = 1_300_000 + + +def q(host: str, graph: str, cypher: str, + timeout: float = 900.0) -> tuple[dict, float]: + """Run one query. Returns (payload, wall seconds). Errors come back HTTP 200.""" + req = urllib.request.Request( + f"{host}/query?graph={graph}", + data=cypher.encode(), + headers={"Content-Type": "application/json"}, + method="POST", + ) + t0 = time.perf_counter() + with urllib.request.urlopen(req, timeout=timeout) as r: + payload = json.loads(r.read()) + return payload, time.perf_counter() - t0 + + +#: Queries known to be pathological are run ONCE, not `repeat` times. The +#: comma-separated join is in the set deliberately: it is measured to show how +#: bad the wrong phrasing is, and at corpus scale it is a genuine cross product +#: over millions of atoms, so repeating it costs tens of minutes for a number +#: whose whole point is that it is large. +SLOW_ONCE = ("JOIN ORDER: acid pairs, comma-separated",) + + +def timed(host: str, graph: str, cypher: str, n: int = 3, + timeout_s: float = 900.0) -> dict: + server_ms: list[float] = [] + wall_ms: list[float] = [] + err = "" + for _ in range(n): + try: + payload, wall = q(host, graph, cypher, timeout_s) + except Exception as exc: # noqa: BLE001 + return {"error": f"{type(exc).__name__}: {exc}"[:160]} + if payload.get("error"): + err = f"{payload['error']}: {payload.get('error_details', '')}" + break + server_ms.append(float(payload.get("time", 0.0))) + wall_ms.append(wall * 1000.0) + if err: + return {"error": err} + return { + "server_ms": round(statistics.median(server_ms), 1), + "wall_ms": round(statistics.median(wall_ms), 1), + } + + +def dir_stats(graph_dir: Path) -> dict: + """Datapart count and byte size of one graph's directory.""" + if not graph_dir.exists(): + return {"exists": False} + parts = [p for p in graph_dir.rglob("*") if p.is_dir() and "datapart" in p.name.lower()] + total = sum(p.stat().st_size for p in graph_dir.rglob("*") if p.is_file()) + return { + "exists": True, + "dataparts": len(parts), + "bytes": total, + "human": f"{total / 1e9:.2f} GB", + "entries": sum(1 for _ in graph_dir.iterdir()), + } + + +#: The queries the UI ships, named as the investigation that issues them. +def query_set(graph: str) -> list[tuple[str, str]]: + # The corrected R2,2(8) dimer query: one connected chain, the ring/catemer + # distinction carried by the precomputed `is_involution` flag, and scoped to + # a single structure because Component is deduplicated corpus-wide. + acid_chain = ( + "MATCH (f1:Fragment)<-[:HAS_FRAGMENT]-(c1:Component)<-[:IN_COMPONENT]-(a1:Atom)" + "-[h:CONTACT]->(a2:Atom)-[:IN_COMPONENT]->(c2:Component)-[:HAS_FRAGMENT]->(f2:Fragment) " + "WHERE f1.fragment_type = 'carboxylic_acid' AND f2.fragment_type = 'carboxylic_acid' " + "AND h.kind = 'hbond' AND h.h_inferred = false AND h.is_involution = true " + "AND a1.element = 'O' AND a2.element = 'O' AND a1.cod_id = a2.cod_id " + "RETURN count(h)" + ) + # Same question, written the way most people would write it: separate + # patterns joined in the WHERE. This pair is the join-order demonstration. + acid_comma = ( + "MATCH (a1:Atom)-[h:CONTACT]->(a2:Atom), " + "(a1)-[:IN_COMPONENT]->(c1:Component)-[:HAS_FRAGMENT]->(f1:Fragment), " + "(a2)-[:IN_COMPONENT]->(c2:Component)-[:HAS_FRAGMENT]->(f2:Fragment) " + "WHERE h.kind = 'hbond' AND h.h_inferred = false " + "AND f1.fragment_type = 'carboxylic_acid' " + "AND f2.fragment_type = 'carboxylic_acid' RETURN count(h)" + ) + acid_chain_1hop = ( + "MATCH (f1:Fragment)<-[:HAS_FRAGMENT]-(c1:Component)<-[:IN_COMPONENT]-(a1:Atom)" + "-[h:CONTACT]->(a2:Atom)-[:IN_COMPONENT]->(c2:Component)-[:HAS_FRAGMENT]->(f2:Fragment) " + "WHERE f1.fragment_type = 'carboxylic_acid' " + "AND f2.fragment_type = 'carboxylic_acid' " + "AND h.kind = 'hbond' AND h.h_inferred = false RETURN count(h)" + ) + return [ + ("count structures", "MATCH (s:Structure) RETURN count(s)"), + ("count atoms", "MATCH (n:Atom) RETURN count(n)"), + ("count contacts", "MATCH ()-[r:CONTACT]->() RETURN count(r)"), + ("spacegroup P2(1)/c", + "MATCH (s:Structure)-[:IN_SPACE_GROUP]->(x:SpaceGroup) " + "WHERE x.hm_symbol = 'P 1 21/c 1' RETURN count(s)"), + ("mean H...A geometry", + "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' AND r.h_inferred = false " + "RETURN avg(r.length)"), + ("dimensionality 2D", + "MATCH (s:Structure) WHERE s.net_dim = 2 RETURN count(s)"), + ("packing subgraph (one structure)", + "MATCH (a:Atom)-[e:CONTACT]->(b:Atom) WHERE a.cod_id = 2229029 " + "RETURN a, labels(a), a.label, e, edgeType(e), b, labels(b), b.label"), + ("fragment census", + "MATCH (c:Component)-[:HAS_FRAGMENT]->(f:Fragment) " + "WHERE f.fragment_type = 'carboxylic_acid' RETURN count(c)"), + ("synthon competition (acid->pyridine)", + "MATCH (f1:Fragment)<-[:HAS_FRAGMENT]-(c1:Component)<-[:IN_COMPONENT]-(a1:Atom)" + "-[h:CONTACT]->(a2:Atom)-[:IN_COMPONENT]->(c2:Component)-[:HAS_FRAGMENT]->(f2:Fragment) " + "WHERE f1.fragment_type = 'carboxylic_acid' " + "AND f2.fragment_type = 'pyridine_nitrogen' AND h.kind = 'hbond' " + "AND h.h_inferred = false AND a1.element = 'O' AND a2.element = 'N' " + "RETURN count(h)"), + ("JOIN ORDER: acid pairs, one chain", acid_chain_1hop), + ("JOIN ORDER: acid pairs, comma-separated", acid_comma), + ("acid dimer R2,2(8)", acid_chain), + ] + + +def main(argv: list[str] | None = None) -> int: + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--host", default="http://localhost:6691") + ap.add_argument("--graph", default="cod_slice_v2") + ap.add_argument("--turing-dir", type=Path, default=Path("turing-data"), + help="graph directory, for the disk and CSD-scale figures") + ap.add_argument("--repeat", type=int, default=3) + ap.add_argument("--merge", action="store_true", + help="also run MERGE_DATAPARTS and re-measure") + ap.add_argument("--out", type=Path, + default=Path("data/manifest/perf.json")) + args = ap.parse_args(argv) + + host, graph = args.host, args.graph + gdir = args.turing_dir / "graphs" / graph + report: dict = {"graph": graph, "host": host} + + print(f"== {graph} ==\n") + before = dir_stats(gdir) + report["disk_before"] = before + if before.get("exists"): + print(f"on disk: {before['human']} in {before['entries']} entries " + f"({before['dataparts']} datapart dirs)\n") + + qs = query_set(graph) + print(f"{'query':<42}{'server ms':>11}{'wall ms':>10}") + print("-" * 63) + results: dict = {} + for name, cypher in qs: + r = timed(host, graph, cypher, + 1 if name in SLOW_ONCE else args.repeat) + results[name] = r + if "error" in r: + print(f"{name:<42}{'ERROR':>11} {r['error'][:40]}") + else: + print(f"{name:<42}{r['server_ms']:>11.1f}{r['wall_ms']:>10.1f}") + report["queries"] = results + + # the join-order headline, stated explicitly + chain = results.get("JOIN ORDER: acid pairs, one chain", {}) + comma = results.get("JOIN ORDER: acid pairs, comma-separated", {}) + if "server_ms" in chain and "server_ms" in comma and chain["server_ms"] > 0: + ratio = comma["server_ms"] / chain["server_ms"] + report["join_order_speedup"] = round(ratio, 1) + print(f"\njoin order: comma-separated is {ratio:.1f}x slower than one " + f"connected chain for the identical question") + + # scale + ns, _ = q(host, graph, "MATCH (s:Structure) RETURN count(s)") + n_struct = ns["data"][0][0][0] if not ns.get("error") else 0 + report["structures"] = n_struct + if before.get("exists") and n_struct: + per = before["bytes"] / n_struct + report["bytes_per_structure"] = round(per) + report["csd_projection_bytes"] = round(per * CSD_STRUCTURES) + print(f"\nscale: {n_struct:,} structures, {per / 1024:.1f} kB per structure " + f"on disk\n -> the CSD at {CSD_STRUCTURES:,} structures " + f"projects to {per * CSD_STRUCTURES / 1e9:.0f} GB") + + if args.merge: + print("\n== MERGE_DATAPARTS ==") + chg, _ = q(host, graph, "CHANGE LIST") + print(f"open changes before merge: {chg.get('data')}") + t0 = time.perf_counter() + m, _ = q(host, graph, "MERGE_DATAPARTS") + print(f"returned in {(time.perf_counter() - t0) * 1000:.1f} ms: " + f"{m.get('error') or 'no error'}") + time.sleep(2) + after = dir_stats(gdir) + report["disk_after_merge"] = after + if after.get("exists"): + print(f"on disk: {before['human']} -> {after['human']} " + f"({before['dataparts']} -> {after['dataparts']} datapart dirs, " + f"{before['entries']} -> {after['entries']} entries)") + post: dict = {} + for name, cypher in qs: + post[name] = timed(host, graph, cypher, + 1 if name in SLOW_ONCE else args.repeat) + report["queries_after_merge"] = post + # correctness, not just speed: the 1.35 run returned CORRUPT reads + changed = [ + n for n in results + if results[n].get("error") != post[n].get("error") + ] + print(f"queries whose outcome changed after merge: {len(changed)}") + + args.out.parent.mkdir(parents=True, exist_ok=True) + args.out.write_text(json.dumps(report, indent=2) + "\n") + print(f"\nwrote {args.out}") + + # Say which promised outputs are missing, and fail. Silently omitting the + # disk figures because --turing-dir pointed nowhere, or the join-order + # ratio because one side timed out, leaves a report that looks complete. + absent = [k for k in ("join_order_speedup", "bytes_per_structure", + "csd_projection_bytes") if k not in report] + errored = [n for n, r in results.items() if "error" in r] + if absent or errored: + print() + if absent: + print("MISSING from the report: " + ", ".join(absent)) + if not report.get("disk_before", {}).get("exists"): + print(f" ({args.turing_dir} does not hold this graph -- " + f"pass --turing-dir)") + if errored: + print("queries that did not complete: " + ", ".join(errored)) + return 1 + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/probe_binary_choice.py b/tests/probe_binary_choice.py deleted file mode 100644 index a566583..0000000 --- a/tests/probe_binary_choice.py +++ /dev/null @@ -1,72 +0,0 @@ -"""Decide which ALREADY-BUILT turingdb binary this demo should run on. - -Two candidates are already on the box (no new build): - A. ~/turingdb-src/build/turing_install/bin/turingdb (source build 3c113474b, - 2026-07-16) - the one proj2 and def_demo drive the visualizer with. - B. .venv/.../turingdb/bin/turingdb (v1.37 wheel, 2026-08-21). - -The demo is being grafted onto the visualizer, so the binary must satisfy BOTH: - 1. the visualizer's data layer -> db.getNodes / db.getNodeEdges / db.getEdges - 2. Q2's deep traversal -> postfix quantifiers -[e]->{1,8} - -Run against each candidate and compare. -""" - -from __future__ import annotations - -import os -import sys - -from turingdb import TuringDB - -HOST = os.environ.get("TURING_HOST", "http://localhost:6691") -GRAPH = os.environ.get("PROBE_GRAPH", "probe_critical") - -CHECKS: list[tuple[str, str]] = [ - # --- visualizer data layer ------------------------------------------- - ("db.getNodes(['Atom'])", "CALL db.getNodes(['Atom'])"), - ("db.getNodes single label", "CALL db.getNodes('Atom')"), - ("db.getNodes with limit", "CALL db.getNodes(['Atom'], 10)"), - ("db.getNodeEdges", "CALL db.getNodeEdges([0], 10, [], [], [], [], false)"), - ("db.getEdges", "CALL db.getEdges(['CONTACT'])"), - ("db.listNodes", "CALL db.listNodes()"), - ("LIST AVAILABLE GRAPHS", "LIST AVAILABLE GRAPHS"), - ("LIST GRAPH", "LIST GRAPH"), - # --- Q2 traversal ------------------------------------------------------ - ("quantifier {1,8}", "MATCH (a:Atom {uid:'A0'})-[e]->{1,8}(b:Atom) RETURN count(b)"), - ("quantifier {8,8}", "MATCH (a:Atom {uid:'A0'})-[e]->{8,8}(b:Atom) RETURN count(b)"), - ("quantifier +", "MATCH (a:Atom {uid:'A0'})-[e]->+(b:Atom) RETURN count(b)"), - # --- versioning -------------------------------------------------------- - ("CALL db.history()", "CALL db.history()"), - ("CALL db.showIndexes()", "CALL db.showIndexes()"), -] - - -def main() -> int: - print(f"--- {HOST} ---") - client = TuringDB(host=HOST) - try: - client.query(f"LOAD GRAPH {GRAPH}") - except Exception: # noqa: BLE001 - pass - client.set_graph(GRAPH) - - for label, cypher in CHECKS: - try: - df = client.query(cypher) - except Exception as exc: # noqa: BLE001 - msg = " ".join(str(exc).split()) - # a dead socket here means the statement killed the server - dead = "Connection refused" in msg or "peer closed" in msg or "incomplete" in msg - tag = "CRASH" if dead else "FAIL " - print(f" {tag} {label:<26} {msg[:95]}") - if dead: - print(" !! server appears to have died - aborting") - return 1 - else: - print(f" PASS {label:<26} rows={len(df)} cols={list(df.columns)[:4]}") - return 0 - - -if __name__ == "__main__": - sys.exit(main()) diff --git a/tests/probe_critical.py b/tests/probe_critical.py deleted file mode 100644 index d6a3ef2..0000000 --- a/tests/probe_critical.py +++ /dev/null @@ -1,250 +0,0 @@ -"""Fast, focused probe of the Cypher constructs the demo cannot do without. - -Separate from probe_dialect.py so it can be pointed at a candidate server build -in a few seconds. Builds a 12-atom hydrogen-bond chain and asks: can this engine -express the four demo queries at all? - -Usage: - TURING_HOST=http://localhost:6691 ./.venv/bin/python tests/probe_critical.py -""" - -from __future__ import annotations - -import os -import sys - -from turingdb import TuringDB, TuringDBException - -HOST = os.environ.get("TURING_HOST", "http://localhost:6691") -GRAPH = os.environ.get("PROBE_GRAPH", "probe_critical") -N_CHAIN = 12 - -# (label, cypher, why it matters) -CRITICAL: list[tuple[str, str, str]] = [ - # --- Q2, the headline "deep multi-hop traversal" claim ----------------- - ("var-len *1..3", - "MATCH (a:Atom {uid:'A0'})-[:CONTACT*1..3]->(b:Atom) RETURN count(b)", - "Q2 as specified"), - ("var-len *1..8", - "MATCH (a:Atom {uid:'A0'})-[:CONTACT*1..8]->(b:Atom) RETURN count(b)", - "Q2 at full depth"), - ("var-len + inline rel prop", - "MATCH (a:Atom {uid:'A0'})-[:CONTACT*1..8 {kind:'hbond'}]->(b:Atom) RETURN count(b)", - "Q2 restricted to hydrogen bonds"), - ("named path + length()", - "MATCH p = (a:Atom {uid:'A0'})-[:CONTACT*1..4]->(b:Atom) RETURN length(p), count(b)", - "Q2 growth curve (reached-per-depth)"), - ("explicit 8-hop chain (fallback)", - "MATCH (a:Atom {uid:'A0'})-[:CONTACT]->(b1:Atom)-[:CONTACT]->(b2:Atom)" - "-[:CONTACT]->(b3:Atom)-[:CONTACT]->(b4:Atom)-[:CONTACT]->(b5:Atom)" - "-[:CONTACT]->(b6:Atom)-[:CONTACT]->(b7:Atom)-[:CONTACT]->(b8:Atom) RETURN count(b8)", - "Q2 fallback if var-len is absent"), - ("shortestPath", - "MATCH (a:Atom {uid:'A0'}), (b:Atom {uid:'A5'}) RETURN shortestPath((a)-[:CONTACT*]->(b))", - "nice-to-have"), - - # --- Q1, the synthon cycle -------------------------------------------- - ("reciprocal cycle pattern", - "MATCH (a1:Atom)-[h1:CONTACT]->(a2:Atom), (a2)-[h2:CONTACT]->(a1) " - "WHERE h1.length < 2.9 AND h2.length < 2.9 RETURN a1.uid, a2.uid, h1.length", - "Q1 acid-dimer ring"), - ("cycle via two MATCH clauses", - "MATCH (a1:Atom)-[h1:CONTACT]->(a2:Atom) " - "MATCH (a2)-[h2:CONTACT]->(a1) RETURN a1.uid, a2.uid", - "Q1 alternative phrasing"), - - # --- Q3, coformer discovery ------------------------------------------- - ("EXISTS { } subquery", - "MATCH (c:Component) WHERE EXISTS { MATCH (c)<-[:CONTAINS_COMPONENT]-(:Structure) } " - "RETURN c.name", - "Q3 'not yet co-crystallised with' clause"), - ("NOT EXISTS { } subquery", - "MATCH (c:Component) WHERE NOT EXISTS { MATCH (c)-[:HAS_FRAGMENT]->(:Fragment) } " - "RETURN c.name", - "Q3 negation"), - ("OPTIONAL MATCH", - "MATCH (s:Structure) OPTIONAL MATCH (s)-[:SUPERSEDES]->(o:Structure) " - "RETURN s.cod_id, o.cod_id", - "Q3 fallback for NOT EXISTS; audit view"), - ("5-hop similarity chain", - "MATCH (api:Component)-[:HAS_FRAGMENT]->(f:Fragment)<-[:HAS_FRAGMENT]-(sim:Component)" - "<-[:CONTAINS_COMPONENT]-(st:Structure)-[:CONTAINS_COMPONENT]->(cand:Component) " - "RETURN cand.name, count(cand)", - "Q3 core traversal"), - - # --- Q4 + reporting ---------------------------------------------------- - ("aggregate + non-aggregate", - "MATCH (f:Fragment)<-[:HAS_FRAGMENT]-(c:Component) RETURN f.fragment_type, count(c)", - "Q4 motif-frequency table"), - ("ORDER BY aggregate", - "MATCH (f:Fragment)<-[:HAS_FRAGMENT]-(c:Component) RETURN f.fragment_type, count(c) AS n ORDER BY n DESC", - "Q4 ranking"), - ("WITH", - "MATCH (a:Atom) WITH count(a) AS n RETURN n", - "query composition"), - ("collect()", - "MATCH (a:Atom) RETURN collect(a.uid)", - "packing-network assembly for the UI"), - ("DISTINCT", - "MATCH (s:Structure) RETURN DISTINCT s.hm_symbol", - "dedup"), - ("count(DISTINCT x)", - "MATCH (s:Structure) RETURN count(DISTINCT s.hm_symbol)", - "Q3 support counting"), - ("UNWIND", - "UNWIND [1,2,3] AS x RETURN x", - "seed injection"), - ("IN list literal", - "MATCH (f:Fragment) WHERE f.fragment_type IN ['carboxylic_acid','hydroxyl'] RETURN f.fragment_type", - "seed injection"), - - # --- general ergonomics ----------------------------------------------- - ("two sequential MATCH clauses", - "MATCH (s:Structure)-[:HAS_SITE]->(a:Atom) MATCH (a)-[e:CONTACT]->(b:Atom) " - "RETURN s.cod_id, count(e)", - "query readability"), - ("IS NOT NULL", - "MATCH (s:Structure) WHERE s.temperature IS NOT NULL RETURN s.cod_id", - "COD has many missing fields"), - ("min/max", - "MATCH ()-[e:CONTACT]->() RETURN min(e.length), max(e.length)", - "reporting"), - ("avg", - "MATCH ()-[e:CONTACT]->() RETURN avg(e.length)", - "reporting"), - ("rel type alternation", - "MATCH (a)-[e:CONTACT|BONDED_TO]->(b) RETURN count(e)", - "convenience"), - ("STARTS WITH", - "MATCH (s:Structure) WHERE s.formula STARTS WITH 'C8' RETURN s.cod_id", - "formula search in Explore screen"), - ("CONTAINS", - "MATCH (s:Structure) WHERE s.formula CONTAINS 'O4' RETURN s.cod_id", - "formula search in Explore screen"), - ("CASE", - "MATCH (s:Structure) RETURN CASE WHEN s.r_factor < 0.04 THEN 'good' ELSE 'ok' END", - "convenience"), - ("id()", - "MATCH (a:Atom) RETURN id(a) LIMIT 2", - "internal id access"), - ("labels()", - "MATCH (n) RETURN labels(n) LIMIT 3", - "UI node typing"), - ("keys()", - "MATCH (s:Structure) RETURN keys(s)", - "UI property panel"), - ("parameters $x", - "MATCH (a:Atom {uid:$seed}) RETURN a.uid", - "safe seed injection"), -] - - -def build(client: TuringDB) -> None: - try: - client.create_graph(GRAPH) - except TuringDBException: - pass - try: - client.query(f"LOAD GRAPH {GRAPH}") - except Exception: # noqa: BLE001 - pass - client.set_graph(GRAPH) - - try: - client.new_change() - client.query("MATCH (n) DETACH DELETE n") - client.query("CHANGE SUBMIT") - client.checkout() - except Exception: # noqa: BLE001 - pass - - nodes = [ - f"CREATE (:Atom {{uid:'A{i}', label:'O{i}', element:'O', occupancy:1.0}})" - for i in range(N_CHAIN) - ] + [ - "CREATE (:Structure {cod_id:1000001, r_factor:0.041, temperature:293.0, " - "formula:'C8 H6 O4', hm_symbol:'P 21/c'})", - "CREATE (:Structure {cod_id:1000002, r_factor:0.038, temperature:100.0, " - "formula:'C8 H6 O4', hm_symbol:'P 21/c'})", - "CREATE (:Component {inchikey:'KDY', formula:'C8 H6 O4', is_solvent:false, " - "name:'terephthalic acid'})", - "CREATE (:Component {inchikey:'XLY', formula:'H2 O', is_solvent:true, name:'water'})", - "CREATE (:Fragment {fragment_type:'carboxylic_acid', smarts:'[CX3](=O)[OX2H1]'})", - "CREATE (:Fragment {fragment_type:'aromatic_ring', smarts:'c1ccccc1'})", - ] - client.new_change() - client.query("\n".join(nodes)) - client.query("COMMIT") # a MATCH cannot see CREATEs from the same change - - edges = [] - for i in range(N_CHAIN - 1): - edges.append( - f"MATCH (x:Atom {{uid:'A{i}'}}), (y:Atom {{uid:'A{i + 1}'}}) " - f"CREATE (x)-[:CONTACT {{kind:'hbond', length:{2.0 + 0.01 * i:.3f}, " - f"angle:165.0, symop:'x,y,z', h_inferred:false}}]->(y)" - ) - edges += [ - "MATCH (x:Atom {uid:'A1'}), (y:Atom {uid:'A0'}) " - "CREATE (x)-[:CONTACT {kind:'hbond', length:2.05, angle:170.0, " - "symop:'-x,-y,-z', h_inferred:false}]->(y)", - "MATCH (x:Atom {uid:'A0'}), (y:Atom {uid:'A1'}) " - "CREATE (x)-[:BONDED_TO {length:1.43, bond_order:1}]->(y)", - "MATCH (s:Structure {cod_id:1000001}), (c:Component {is_solvent:false}) " - "CREATE (s)-[:CONTAINS_COMPONENT {n_copies:2, role:'principal'}]->(c)", - "MATCH (s:Structure {cod_id:1000001}), (c:Component {is_solvent:true}) " - "CREATE (s)-[:CONTAINS_COMPONENT {n_copies:1, role:'solvent'}]->(c)", - "MATCH (s:Structure {cod_id:1000002}), (c:Component {is_solvent:false}) " - "CREATE (s)-[:CONTAINS_COMPONENT {n_copies:1, role:'principal'}]->(c)", - "MATCH (c:Component {is_solvent:false}), (f:Fragment {fragment_type:'carboxylic_acid'}) " - "CREATE (c)-[:HAS_FRAGMENT {atom_uids:'A0,A1'}]->(f)", - "MATCH (c:Component {is_solvent:true}), (f:Fragment {fragment_type:'aromatic_ring'}) " - "CREATE (c)-[:HAS_FRAGMENT {atom_uids:'A4,A5'}]->(f)", - "MATCH (s:Structure {cod_id:1000001}), (a:Atom {uid:'A0'}) CREATE (s)-[:HAS_SITE]->(a)", - "MATCH (s:Structure {cod_id:1000001}), (a:Atom {uid:'A1'}) CREATE (s)-[:HAS_SITE]->(a)", - "MATCH (x:Structure {cod_id:1000002}), (y:Structure {cod_id:1000001}) " - "CREATE (x)-[:SUPERSEDES {reason:'redetermined at 100K', inferred:true}]->(y)", - ] - for stmt in edges: - client.query(stmt) - client.query("CHANGE SUBMIT") - client.checkout() - - -def main() -> int: - print(f"connecting to {HOST}") - client = TuringDB(host=HOST) - print("building fixture ...") - build(client) - n = client.query("MATCH (n) RETURN count(n)").iloc[0, 0] - e = client.query("MATCH ()-[r]->() RETURN count(r)").iloc[0, 0] - print(f"fixture: {n} nodes, {e} edges\n") - - width = max(len(name) for name, _, _ in CRITICAL) - npass = 0 - fails: list[tuple[str, str, str]] = [] - for name, cypher, why in CRITICAL: - try: - df = client.query(cypher) - except Exception as exc: # noqa: BLE001 - msg = str(exc).replace("\n", " ") - # keep just the engine's reason, drop the echoed query block - for marker in ("-------* ", ": "): - if marker in msg: - pass - short = msg[:150] - print(f" FAIL {name:<{width}} {short}") - fails.append((name, why, short)) - else: - npass += 1 - print(f" PASS {name:<{width}} rows={len(df)} cols={list(df.columns)}") - - print(f"\n{npass}/{len(CRITICAL)} critical constructs supported") - if fails: - print("\nMISSING (impact):") - for name, why, msg in fails: - print(f" - {name} -> {why}") - return 0 - - -if __name__ == "__main__": - sys.exit(main()) diff --git a/tests/probe_dialect.py b/tests/probe_dialect.py deleted file mode 100644 index aae584a..0000000 --- a/tests/probe_dialect.py +++ /dev/null @@ -1,532 +0,0 @@ -"""Probe the TuringDB Cypher dialect against a live instance. - -Run BEFORE writing any real query code. Builds a tiny graph shaped like the -crystallography schema, then tries every Cypher construct the demo needs and -records PASS/FAIL with the server's own error text. - -Usage: - ./.venv/bin/python tests/probe_dialect.py # human-readable report - ./.venv/bin/python tests/probe_dialect.py --json # machine-readable - -Nothing here is load-bearing for the demo; it exists so that dialect gaps show -up on day one rather than after 500 lines of Cypher have been written. -""" - -from __future__ import annotations - -import argparse -import json -import os -import sys -import traceback -from dataclasses import dataclass, field -from typing import Any, Callable - -from turingdb import TuringDB, TuringDBException - -HOST = os.environ.get("TURING_HOST", "http://localhost:6690") -GRAPH = os.environ.get("PROBE_GRAPH", "probe_dialect") - - -# -------------------------------------------------------------------------- -# result bookkeeping -# -------------------------------------------------------------------------- - - -@dataclass -class Result: - name: str - group: str - cypher: str - ok: bool - detail: str = "" - rows: int | None = None - columns: list[str] = field(default_factory=list) - sample: str = "" - - -RESULTS: list[Result] = [] - - -def _short(exc: BaseException) -> str: - """Collapse a TuringDB error to a single readable line.""" - text = str(exc).strip().replace("\n", " ") - return text[:400] if text else exc.__class__.__name__ - - -def probe(name: str, group: str, cypher: str, *, client: TuringDB) -> Result: - """Run one read query and record whether the dialect accepted it.""" - try: - df = client.query(cypher) - except (TuringDBException, Exception) as exc: # noqa: BLE001 - want everything - res = Result(name, group, cypher, ok=False, detail=_short(exc)) - else: - cols = [str(c) for c in getattr(df, "columns", [])] - sample = "" - try: - if len(df) > 0: - sample = df.head(2).to_string(index=False).replace("\n", " | ")[:220] - except Exception: # noqa: BLE001 - sample = "" - res = Result(name, group, cypher, ok=True, rows=len(df), columns=cols, sample=sample) - RESULTS.append(res) - return res - - -def probe_call(name: str, group: str, fn: Callable[[], Any], *, label: str = "") -> Result: - """Run one SDK call (not raw Cypher) and record whether it worked.""" - try: - out = fn() - except Exception as exc: # noqa: BLE001 - res = Result(name, group, label or f" {name}", ok=False, detail=_short(exc)) - else: - text = repr(out) - if len(text) > 300: - text = text[:300] + "..." - res = Result(name, group, label or f" {name}", ok=True, sample=text) - RESULTS.append(res) - return res - - -# -------------------------------------------------------------------------- -# fixture graph -# -------------------------------------------------------------------------- - -# A chain long enough to exercise 8-hop traversal, plus a two-way cycle that -# mimics the carboxylic-acid dimer of Q1, plus one structure/component/fragment -# spine so the join-heavy queries have something to bind to. -# -# Property names are chosen to match the real schema EXACTLY, because names -# like `length`, `order`, `angle`, `count` and `position` collide with Cypher -# keywords/functions in some dialects and we need to find that out now. - -N_CHAIN = 12 - - -def build_fixture(client: TuringDB) -> None: - """(Re)build the probe graph from scratch. Idempotent.""" - try: - client.create_graph(GRAPH) - print(f" created graph {GRAPH}") - except TuringDBException: - print(f" graph {GRAPH} already exists") - try: - client.load_graph(GRAPH) - except TuringDBException: - pass - try: - client.query(f"LOAD GRAPH {GRAPH}") - except Exception: # noqa: BLE001 - pass - client.set_graph(GRAPH) - - # wipe - try: - client.new_change() - client.query("MATCH (n) DETACH DELETE n") - client.query("CHANGE SUBMIT") - client.checkout() - except Exception as exc: # noqa: BLE001 - print(f" wipe skipped: {_short(exc)}") - - # --- nodes ------------------------------------------------------------- - atom_stmts = [] - for i in range(N_CHAIN): - atom_stmts.append( - f"CREATE (:Atom {{uid:'A{i}', label:'O{i}', element:'O', " - f"fract_x:{0.1 * i:.3f}, fract_y:0.2, fract_z:0.3, " - f"occupancy:1.0, u_iso:0.05}})" - ) - node_cypher = "\n".join( - atom_stmts - + [ - "CREATE (:Structure {cod_id:1000001, a:5.1, b:6.2, c:7.3, alpha:90.0, " - "beta:99.5, gamma:90.0, volume:230.4, Z:4, Z_prime:1.0, r_factor:0.041, " - "temperature:293.0, pressure:101.325, formula:'C8 H6 O4', " - "hm_symbol:'P 21/c', revision_kind:'redetermination'})", - "CREATE (:Structure {cod_id:1000002, a:5.2, b:6.3, c:7.4, alpha:90.0, " - "beta:99.9, gamma:90.0, volume:231.9, Z:4, Z_prime:1.0, r_factor:0.038, " - "temperature:100.0, pressure:101.325, formula:'C8 H6 O4', " - "hm_symbol:'P 21/c', revision_kind:'correction'})", - "CREATE (:SpaceGroup {hm_symbol:'P 21/c', number:14, crystal_system:'monoclinic'})", - "CREATE (:Element {symbol:'O', atomic_number:8})", - "CREATE (:Component {inchikey:'KDYFGRWQOYBRFD-UHFFFAOYSA-N', formula:'C8 H6 O4', " - "charge:0, is_solvent:false, name:'terephthalic acid'})", - "CREATE (:Component {inchikey:'XLYOFNOQVPJJNP-UHFFFAOYSA-N', formula:'H2 O', " - "charge:0, is_solvent:true, name:'water'})", - "CREATE (:Fragment {type:'carboxylic_acid', smarts:'[CX3](=O)[OX2H1]'})", - "CREATE (:Fragment {type:'aromatic_ring', smarts:'c1ccccc1'})", - # NB: `volume` is deliberately NOT reused here. TuringDB types - # properties GLOBALLY, not per-label, so Publication.volume (String) - # collides with Structure.volume (Double, the cell volume) and the - # whole CREATE is rejected. Hence `journal_volume`. - "CREATE (:Publication {doi:'10.1039/x', year:2015, title:'A paper', " - "journal_volume:'17', pages:'1-9'})", - "CREATE (:Author {name:'A. Crystallographer', name_key:'crystallographer_a'})", - "CREATE (:Journal {name:'CrystEngComm', issn:'1466-8033'})", - "CREATE (:Deposition {deposition_id:'D1', deposited_date:'2015-03-04', " - "depositor:'dep'})", - ] - ) - client.new_change() - client.query(node_cypher) - # A MATCH inside a change cannot see nodes CREATEd in the same change, so - # the node writes must be committed before the edge writes can bind them. - client.query("COMMIT") - - # --- edges ------------------------------------------------------------- - edge_stmts = [] - # hbond chain A0 -> A1 -> ... for deep traversal - for i in range(N_CHAIN - 1): - edge_stmts.append( - f"MATCH (x:Atom {{uid:'A{i}'}}), (y:Atom {{uid:'A{i + 1}'}}) " - f"CREATE (x)-[:CONTACT {{kind:'hbond', length:{2.0 + 0.01 * i:.3f}, " - f"angle:165.0, symop:'x,y,z', h_inferred:false}}]->(y)" - ) - # reciprocal pair => the Q1 acid-dimer cycle - edge_stmts.append( - "MATCH (x:Atom {uid:'A0'}), (y:Atom {uid:'A1'}) " - "CREATE (y)-[:CONTACT {kind:'hbond', length:2.05, angle:170.0, " - "symop:'-x,-y,-z', h_inferred:false}]->(x)" - ) - edge_stmts.append( - "MATCH (x:Atom {uid:'A2'}), (y:Atom {uid:'A3'}) " - "CREATE (x)-[:CONTACT {kind:'halogen', length:3.30, angle:158.0, " - "symop:'x,1/2-y,z', h_inferred:false}]->(y)" - ) - edge_stmts.append( - "MATCH (x:Atom {uid:'A4'}), (y:Atom {uid:'A5'}) " - "CREATE (x)-[:CONTACT {kind:'pi_stack', length:3.62, angle:8.0, " - "symop:'x,y,1+z', h_inferred:true}]->(y)" - ) - edge_stmts.append( - "MATCH (x:Atom {uid:'A0'}), (y:Atom {uid:'A1'}) " - "CREATE (x)-[:BONDED_TO {length:1.43, order:1}]->(y)" - ) - edge_stmts += [ - "MATCH (s:Structure {cod_id:1000001}), (g:SpaceGroup {number:14}) " - "CREATE (s)-[:IN_SPACE_GROUP]->(g)", - "MATCH (s:Structure {cod_id:1000001}), (e:Element {symbol:'O'}) " - "CREATE (s)-[:CONTAINS_ELEMENT]->(e)", - "MATCH (s:Structure {cod_id:1000001}), (c:Component {is_solvent:false}) " - "CREATE (s)-[:CONTAINS_COMPONENT {count:2, role:'principal'}]->(c)", - "MATCH (s:Structure {cod_id:1000001}), (c:Component {is_solvent:true}) " - "CREATE (s)-[:CONTAINS_COMPONENT {count:1, role:'solvent'}]->(c)", - "MATCH (s:Structure {cod_id:1000002}), (c:Component {is_solvent:false}) " - "CREATE (s)-[:CONTAINS_COMPONENT {count:1, role:'principal'}]->(c)", - "MATCH (c:Component {is_solvent:false}), (f:Fragment {type:'carboxylic_acid'}) " - "CREATE (c)-[:HAS_FRAGMENT {atom_uids:'A0,A1'}]->(f)", - "MATCH (c:Component {is_solvent:false}), (f:Fragment {type:'aromatic_ring'}) " - "CREATE (c)-[:HAS_FRAGMENT {atom_uids:'A4,A5'}]->(f)", - "MATCH (f:Fragment {type:'carboxylic_acid'}), (a:Atom {uid:'A0'}) " - "CREATE (f)-[:INCLUDES]->(a)", - "MATCH (f:Fragment {type:'carboxylic_acid'}), (a:Atom {uid:'A1'}) " - "CREATE (f)-[:INCLUDES]->(a)", - "MATCH (s:Structure {cod_id:1000001}), (a:Atom {uid:'A0'}) " - "CREATE (s)-[:HAS_SITE]->(a)", - "MATCH (s:Structure {cod_id:1000001}), (a:Atom {uid:'A1'}) " - "CREATE (s)-[:HAS_SITE]->(a)", - "MATCH (s:Structure {cod_id:1000001}), (p:Publication {year:2015}) " - "CREATE (s)-[:PUBLISHED_IN]->(p)", - "MATCH (p:Publication {year:2015}), (j:Journal {issn:'1466-8033'}) " - "CREATE (p)-[:IN_JOURNAL]->(j)", - "MATCH (p:Publication {year:2015}), (au:Author {name_key:'crystallographer_a'}) " - "CREATE (p)-[:AUTHORED_BY {position:1}]->(au)", - "MATCH (s:Structure {cod_id:1000001}), (d:Deposition {deposition_id:'D1'}) " - "CREATE (s)-[:FROM_DEPOSITION]->(d)", - "MATCH (a:Structure {cod_id:1000002}), (b:Structure {cod_id:1000001}) " - "CREATE (a)-[:SUPERSEDES {reason:'redetermined at 100K', inferred:true}]->(b)", - ] - for stmt in edge_stmts: - client.query(stmt) - client.query("CHANGE SUBMIT") - client.checkout() - - n = client.query("MATCH (n) RETURN count(n)") - print(f" fixture built: {n.iloc[0, 0]} nodes") - - -# -------------------------------------------------------------------------- -# the probes -# -------------------------------------------------------------------------- - - -def run_probes(client: TuringDB) -> None: - p = lambda name, group, cypher: probe(name, group, cypher, client=client) # noqa: E731 - - # ---- basics ---------------------------------------------------------- - g = "basics" - p("MATCH label + props", g, "MATCH (s:Structure) RETURN s.cod_id, s.r_factor") - p("inline node prop filter", g, "MATCH (s:Structure {cod_id:1000001}) RETURN s.formula") - p("WHERE on node prop", g, "MATCH (s:Structure) WHERE s.r_factor < 0.04 RETURN s.cod_id") - p("WHERE label predicate", g, "MATCH (n) WHERE n:Atom RETURN count(n)") - p("undirected match", g, "MATCH (a:Atom)-[e:CONTACT]-(b:Atom) RETURN count(e)") - p("return node var", g, "MATCH (s:Structure) RETURN s") - p("return edge var", g, "MATCH (a)-[e:CONTACT]->(b) RETURN a, e, b") - p("IS NOT NULL", g, "MATCH (s:Structure) WHERE s.temperature IS NOT NULL RETURN s.cod_id") - p("boolean prop filter", g, "MATCH (c:Component) WHERE c.is_solvent = false RETURN c.name") - p("OR chain", g, "MATCH (s:Structure) WHERE s.cod_id = 1000001 OR s.cod_id = 1000002 RETURN s.cod_id") - - # ---- property names that may collide with keywords ------------------- - g = "reserved-property-names" - p("rel prop .length", g, "MATCH ()-[e:CONTACT]->() RETURN e.length") - p("rel prop .order", g, "MATCH ()-[e:BONDED_TO]->() RETURN e.order") - p("rel prop .angle", g, "MATCH ()-[e:CONTACT]->() RETURN e.angle") - p("rel prop .count", g, "MATCH ()-[e:CONTAINS_COMPONENT]->() RETURN e.count") - p("rel prop .position", g, "MATCH ()-[e:AUTHORED_BY]->() RETURN e.position") - p("rel prop .role", g, "MATCH ()-[e:CONTAINS_COMPONENT]->() RETURN e.role") - p("rel prop .symop", g, "MATCH ()-[e:CONTACT]->() RETURN e.symop") - p("node prop .number", g, "MATCH (g:SpaceGroup) RETURN g.number") - p("node prop .volume", g, "MATCH (s:Structure) RETURN s.volume") - p("node prop .charge", g, "MATCH (c:Component) RETURN c.charge") - p("node prop .type", g, "MATCH (f:Fragment) RETURN f.type") - p("node prop .name", g, "MATCH (j:Journal) RETURN j.name") - p("node prop .year", g, "MATCH (pb:Publication) RETURN pb.year") - p("node prop .pressure", g, "MATCH (s:Structure) RETURN s.pressure") - p("node prop .Z / .Z_prime", g, "MATCH (s:Structure) RETURN s.Z, s.Z_prime") - p("node prop .a/.b/.c cell", g, "MATCH (s:Structure) RETURN s.a, s.b, s.c") - p("node prop .alpha/.beta/.gamma", g, "MATCH (s:Structure) RETURN s.alpha, s.beta, s.gamma") - p("node prop .element", g, "MATCH (a:Atom) RETURN a.element") - p("node prop .occupancy", g, "MATCH (a:Atom) RETURN a.occupancy") - - # ---- relationship predicates ---------------------------------------- - g = "rel-predicates" - p("inline rel prop filter", g, "MATCH (a)-[e:CONTACT {kind:'hbond'}]->(b) RETURN count(e)") - p("WHERE on rel prop", g, "MATCH (a)-[e:CONTACT]->(b) WHERE e.kind = 'hbond' RETURN count(e)") - p("WHERE rel numeric", g, "MATCH (a)-[e:CONTACT]->(b) WHERE e.length < 2.05 RETURN e.length") - p("WHERE edge label", g, "MATCH (a)-[e]->(b) WHERE e:CONTACT RETURN count(e)") - p("rel type alternation", g, "MATCH (a)-[e:CONTACT|BONDED_TO]->(b) RETURN count(e)") - p("two inline rel props", g, "MATCH (a)-[e:CONTACT {kind:'hbond', h_inferred:false}]->(b) RETURN count(e)") - - # ---- variable length ------------------------------------------------- - g = "variable-length" - p("var-len *1..3", g, "MATCH (a:Atom {uid:'A0'})-[:CONTACT*1..3]->(b:Atom) RETURN count(b)") - p("var-len *1..8", g, "MATCH (a:Atom {uid:'A0'})-[:CONTACT*1..8]->(b:Atom) RETURN count(b)") - p("var-len + rel prop inline", g, - "MATCH (a:Atom {uid:'A0'})-[:CONTACT*1..8 {kind:'hbond'}]->(b:Atom) RETURN count(b)") - p("var-len untyped", g, "MATCH (a:Atom {uid:'A0'})-[*1..2]->(b) RETURN count(b)") - p("var-len exact *3", g, "MATCH (a:Atom {uid:'A0'})-[:CONTACT*3]->(b:Atom) RETURN count(b)") - p("var-len *..8 open lower", g, "MATCH (a:Atom {uid:'A0'})-[:CONTACT*..8]->(b:Atom) RETURN count(b)") - p("named path + length()", g, - "MATCH p = (a:Atom {uid:'A0'})-[:CONTACT*1..4]->(b:Atom) RETURN length(p), count(b)") - p("named path return p", g, "MATCH p = (a:Atom {uid:'A0'})-[:CONTACT*1..2]->(b:Atom) RETURN p") - p("var-len WHERE on rel var", g, - "MATCH (a:Atom {uid:'A0'})-[e:CONTACT*1..4]->(b:Atom) WHERE e.kind = 'hbond' RETURN count(b)") - p("group-by depth", g, - "MATCH p = (a:Atom {uid:'A0'})-[:CONTACT*1..8]->(b:Atom) RETURN length(p), count(b)") - p("shortestPath", g, - "MATCH (a:Atom {uid:'A0'}), (b:Atom {uid:'A5'}) RETURN shortestPath((a)-[:CONTACT*]->(b))") - - # ---- multi-pattern / multi-clause ------------------------------------ - g = "multi-pattern" - p("comma-joined patterns", g, - "MATCH (f:Fragment {type:'carboxylic_acid'})<-[:HAS_FRAGMENT]-(c:Component), " - "(c)<-[:CONTAINS_COMPONENT]-(s:Structure) RETURN s.cod_id, f.type") - p("two sequential MATCH clauses", g, - "MATCH (s:Structure)-[:HAS_SITE]->(a:Atom) " - "MATCH (a)-[e:CONTACT]->(b:Atom) RETURN s.cod_id, count(e)") - p("three sequential MATCH clauses", g, - "MATCH (s:Structure)-[:CONTAINS_COMPONENT]->(c:Component) " - "MATCH (c)-[:HAS_FRAGMENT]->(f:Fragment) " - "MATCH (s)-[:HAS_SITE]->(a:Atom) RETURN s.cod_id, f.type, count(a)") - p("4-hop chain one pattern", g, - "MATCH (s:Structure)-[:CONTAINS_COMPONENT]->(c:Component)-[:HAS_FRAGMENT]->" - "(f:Fragment)<-[:HAS_FRAGMENT]-(c2:Component)<-[:CONTAINS_COMPONENT]-(s2:Structure) " - "RETURN s.cod_id, s2.cod_id, f.type") - p("cycle pattern (Q1 shape)", g, - "MATCH (a1:Atom)-[h1:CONTACT]->(a2:Atom), (a2)-[h2:CONTACT]->(a1) " - "WHERE h1.length < 2.9 AND h2.length < 2.9 RETURN a1.uid, a2.uid, h1.length, h2.length") - p("self-join with inequality", g, - "MATCH (c1:Component)-[:HAS_FRAGMENT]->(f:Fragment)<-[:HAS_FRAGMENT]-(c2:Component) " - "WHERE c1.inchikey <> c2.inchikey RETURN c1.name, c2.name") - - # ---- aggregation / projection --------------------------------------- - g = "aggregation" - p("count(*)", g, "MATCH (a:Atom) RETURN count(*)") - p("count(n)", g, "MATCH (a:Atom) RETURN count(a)") - p("aggregate + non-aggregate", g, - "MATCH (f:Fragment)<-[:HAS_FRAGMENT]-(c:Component) RETURN f.type, count(c)") - p("count with alias", g, "MATCH (a:Atom) RETURN count(a) AS n") - p("collect()", g, "MATCH (a:Atom) RETURN collect(a.uid)") - p("DISTINCT", g, "MATCH (s:Structure) RETURN DISTINCT s.hm_symbol") - p("count(DISTINCT x)", g, "MATCH (s:Structure) RETURN count(DISTINCT s.hm_symbol)") - p("min/max/avg/sum", g, "MATCH (a)-[e:CONTACT]->(b) RETURN min(e.length), max(e.length)") - p("avg", g, "MATCH (a)-[e:CONTACT]->(b) RETURN avg(e.length)") - p("sum", g, "MATCH (a)-[e:CONTACT]->(b) RETURN sum(e.length)") - p("ORDER BY property", g, "MATCH (a)-[e:CONTACT]->(b) RETURN e.length ORDER BY e.length") - p("ORDER BY alias", g, "MATCH (a:Atom) RETURN count(a) AS n ORDER BY n") - p("ORDER BY agg no alias", g, - "MATCH (f:Fragment)<-[:HAS_FRAGMENT]-(c:Component) RETURN f.type, count(c) ORDER BY count(c) DESC") - p("LIMIT", g, "MATCH (a:Atom) RETURN a.uid LIMIT 3") - p("SKIP + LIMIT", g, "MATCH (a:Atom) RETURN a.uid ORDER BY a.uid SKIP 2 LIMIT 3") - - # ---- constructs the spec flagged as uncertain ------------------------ - g = "uncertain-constructs" - p("OPTIONAL MATCH", g, - "MATCH (s:Structure) OPTIONAL MATCH (s)-[:SUPERSEDES]->(old:Structure) " - "RETURN s.cod_id, old.cod_id") - p("EXISTS { } subquery", g, - "MATCH (c:Component) WHERE EXISTS { MATCH (c)<-[:CONTAINS_COMPONENT]-(:Structure) } " - "RETURN c.name") - p("NOT EXISTS { } subquery", g, - "MATCH (c:Component) WHERE NOT EXISTS { MATCH (c)-[:HAS_FRAGMENT]->(:Fragment) } " - "RETURN c.name") - p("negated path predicate", g, - "MATCH (c:Component) WHERE NOT (c)-[:HAS_FRAGMENT]->() RETURN c.name") - p("WITH", g, "MATCH (a:Atom) WITH count(a) AS n RETURN n") - p("WITH + filter", g, - "MATCH (f:Fragment)<-[:HAS_FRAGMENT]-(c:Component) WITH f, count(c) AS n " - "WHERE n > 0 RETURN f.type, n") - p("UNWIND", g, "UNWIND [1,2,3] AS x RETURN x") - p("IN list literal", g, - "MATCH (f:Fragment) WHERE f.type IN ['carboxylic_acid','hydroxyl'] RETURN f.type") - p("keys(n)", g, "MATCH (s:Structure) RETURN keys(s)") - p("labels(n)", g, "MATCH (n) RETURN labels(n) LIMIT 3") - p("labels() in WHERE", g, "MATCH (n) WHERE labels(n) = 'Atom' RETURN count(n)") - p("id() function", g, "MATCH (a:Atom) RETURN id(a) LIMIT 2") - p("CASE expression", g, - "MATCH (s:Structure) RETURN CASE WHEN s.r_factor < 0.04 THEN 'good' ELSE 'ok' END") - p("toInteger", g, "MATCH (s:Structure) WHERE s.cod_id > toInteger('1000000') RETURN s.cod_id") - p("toFloat", g, "MATCH (s:Structure) RETURN s.volume * toFloat('2.0')") - p("arithmetic in RETURN", g, "MATCH (s:Structure) RETURN s.a * s.b * s.c") - p("string equality on prop", g, "MATCH (g2:SpaceGroup) WHERE g2.hm_symbol = 'P 21/c' RETURN g2.number") - p("STARTS WITH", g, "MATCH (s:Structure) WHERE s.formula STARTS WITH 'C8' RETURN s.cod_id") - p("CONTAINS", g, "MATCH (s:Structure) WHERE s.formula CONTAINS 'O4' RETURN s.cod_id") - p("parameter $x", g, "MATCH (a:Atom {uid:$seed}) RETURN a.uid") - p("internal-id seed WHERE a = n", g, "MATCH (a)-->(b) WHERE a = 0 RETURN count(b)") - - # ---- introspection --------------------------------------------------- - g = "introspection" - p("CALL db.labels()", g, "CALL db.labels()") - p("CALL db.edgeTypes()", g, "CALL db.edgeTypes()") - p("CALL db.propertyTypes()", g, "CALL db.propertyTypes()") - p("CALL db.history()", g, "CALL db.history()") - p("CALL db.getNodes", g, "CALL db.getNodes(['Atom'])") - p("LIST GRAPH", g, "LIST GRAPH") - p("LIST AVAILABLE GRAPHS", g, "LIST AVAILABLE GRAPHS") - - -def run_versioning_probes(client: TuringDB) -> None: - """Probe the git-style versioning surface: commits, branches, tags, checkout.""" - g = "versioning-sdk" - - for meth in [ - "new_change", "commit", "checkout", "history", "set_commit", "set_graph", - "create_graph", "load_graph", "list_available_graphs", "list_loaded_graphs", - "create_branch", "branches", "list_branches", "switch_branch", "tag", "tags", - "create_tag", "list_tags", "diff", "merge", "reset", "status", - ]: - has = hasattr(client, meth) - RESULTS.append( - Result(f"SDK has .{meth}()", "versioning-sdk-surface", f"hasattr(client, '{meth}')", - ok=has, detail="" if has else "absent") - ) - - probe_call("client.history()", g, lambda: client.history()) - probe_call("client.list_available_graphs()", g, lambda: client.list_available_graphs()) - probe_call("client.list_loaded_graphs()", g, lambda: client.list_loaded_graphs()) - - g = "versioning-cypher" - pv = lambda name, cypher: probe(name, g, cypher, client=client) # noqa: E731 - pv("CALL db.history()", "CALL db.history()") - pv("LIST BRANCH", "LIST BRANCH") - pv("LIST BRANCHES", "LIST BRANCHES") - pv("SHOW BRANCHES", "SHOW BRANCHES") - pv("LIST TAG", "LIST TAG") - pv("LIST TAGS", "LIST TAGS") - pv("CREATE BRANCH probe_branch", "CREATE BRANCH probe_branch") - pv("BRANCH probe_branch2", "BRANCH probe_branch2") - pv("CREATE TAG probe_tag", "CREATE TAG probe_tag") - pv("TAG probe_tag2", "TAG probe_tag2") - pv("CHECKOUT main", "CHECKOUT main") - - # time travel: grab a real commit hash from history, then try to read at it - try: - hist = client.query("CALL db.history()") - first = str(hist.iloc[0, 0]).replace("(HEAD)", "").strip() - RESULTS.append(Result("history HEAD hash", "versioning-cypher", - "CALL db.history() -> iloc[0,0]", ok=True, sample=first)) - pv(f"LOAD COMMIT ''", f"LOAD COMMIT '{first}'") - probe_call("client.set_commit(hash)", "versioning-sdk", - lambda: client.set_commit(first), label=f"client.set_commit('{first}')") - probe_call("client.checkout() back to HEAD", "versioning-sdk", lambda: client.checkout()) - except Exception as exc: # noqa: BLE001 - RESULTS.append(Result("history HEAD hash", "versioning-cypher", - "CALL db.history()", ok=False, detail=_short(exc))) - - -# -------------------------------------------------------------------------- -# reporting -# -------------------------------------------------------------------------- - - -def report() -> None: - groups: dict[str, list[Result]] = {} - for r in RESULTS: - groups.setdefault(r.group, []).append(r) - - n_ok = sum(1 for r in RESULTS if r.ok) - print("\n" + "=" * 78) - print(f"TuringDB dialect probe — {n_ok}/{len(RESULTS)} constructs accepted") - print(f"host={HOST} graph={GRAPH}") - print("=" * 78) - - for name, rs in groups.items(): - ok = sum(1 for r in rs if r.ok) - print(f"\n### {name} ({ok}/{len(rs)})") - for r in rs: - mark = "PASS" if r.ok else "FAIL" - extra = "" - if r.ok and r.rows is not None: - extra = f" rows={r.rows}" - if r.columns: - extra += f" cols={r.columns}" - print(f" [{mark}] {r.name}{extra}") - if not r.ok: - print(f" cypher: {r.cypher[:150]}") - print(f" error : {r.detail}") - elif r.sample: - print(f" -> {r.sample[:180]}") - - print("\n" + "-" * 78) - print("FAILURES (each of these forces a query rewrite):") - for r in RESULTS: - if not r.ok: - print(f" - [{r.group}] {r.name}: {r.detail[:160]}") - - -def main() -> int: - ap = argparse.ArgumentParser() - ap.add_argument("--json", action="store_true", help="emit JSON instead of a report") - ap.add_argument("--no-build", action="store_true", help="skip fixture rebuild") - args = ap.parse_args() - - print(f"connecting to {HOST} ...") - client = TuringDB(host=HOST) - - if not args.no_build: - print("building fixture graph ...") - build_fixture(client) - else: - try: - client.load_graph(GRAPH) - except Exception: # noqa: BLE001 - pass - client.set_graph(GRAPH) - - print("running probes ...") - run_probes(client) - try: - run_versioning_probes(client) - except Exception: # noqa: BLE001 - traceback.print_exc() - - if args.json: - print(json.dumps([r.__dict__ for r in RESULTS], indent=2)) - else: - report() - return 0 - - -if __name__ == "__main__": - sys.exit(main()) diff --git a/tests/probe_followup.py b/tests/probe_followup.py deleted file mode 100644 index a05bf39..0000000 --- a/tests/probe_followup.py +++ /dev/null @@ -1,156 +0,0 @@ -"""Follow-up probe: nail down the exact workable forms after probe_turingdb_syntax. - -Answers the questions that decide the architecture: - 1. Does an UNTYPED quantified traversal work when seeded? (-> Q2 is possible) - 2. Do IS NULL / IS NOT NULL work when the property exists on some nodes? - (COD is full of missing fields, so this drives the ingest null strategy.) - 3. Do property indexes work, and do they speed up seeded lookups? - 4. Does LOAD JSONL accept the documented APOC shape? -""" - -from __future__ import annotations - -import json -import os -import subprocess -import sys -import time - -from turingdb import TuringDB, TuringDBException - -HOST = os.environ.get("TURING_HOST", "http://localhost:6691") -GRAPH = os.environ.get("PROBE_GRAPH", "probe_critical") -TURING_DIR = os.path.expanduser(os.environ.get("TURING_DIR", "~/ccdc_demo/turing-data-137")) - - -def show(client: TuringDB, label: str, cypher: str) -> None: - try: - df = client.query(cypher) - except Exception as exc: # noqa: BLE001 - print(f" FAIL {label}\n {' '.join(str(exc).split())[:200]}") - else: - val = "" - try: - if len(df) and len(df.columns) == 1: - val = f" = {df.iloc[0, 0]}" - elif len(df): - val = f" first={list(df.iloc[0])[:5]}" - except Exception: # noqa: BLE001 - pass - print(f" PASS {label} rows={len(df)}{val}") - - -def main() -> int: - client = TuringDB(host=HOST) - try: - client.query(f"LOAD GRAPH {GRAPH}") - except Exception: # noqa: BLE001 - pass - client.set_graph(GRAPH) - - print("=== 1. UNTYPED quantified traversal (the Q2 candidate) ===") - show(client, "seeded untyped {1,8}", - "MATCH (a:Atom {uid:'A0'})-[e]->{1,8}(b:Atom) RETURN count(b)") - show(client, "seeded untyped {1,1}", - "MATCH (a:Atom {uid:'A0'})-[e]->{1,1}(b:Atom) RETURN count(b)") - show(client, "seeded untyped {1,4}", - "MATCH (a:Atom {uid:'A0'})-[e]->{1,4}(b:Atom) RETURN count(b)") - show(client, "seeded untyped {8,8} (exact depth)", - "MATCH (a:Atom {uid:'A0'})-[e]->{8,8}(b:Atom) RETURN count(b)") - show(client, "seeded untyped {1,12}", - "MATCH (a:Atom {uid:'A0'})-[e]->{1,12}(b:Atom) RETURN count(b)") - show(client, "untyped {1,8} returning endpoints", - "MATCH (a:Atom {uid:'A0'})-[e]->{1,8}(b:Atom) RETURN b.uid LIMIT 10") - show(client, "untyped one-or-more +", - "MATCH (a:Atom {uid:'A0'})-[e]->+(b:Atom) RETURN count(b)") - show(client, "untyped, seed by internal id", - "MATCH (a)-[e]->{1,8}(b:Atom) WHERE a = 0 RETURN count(b)") - show(client, "untyped {1,8} + WHERE on endpoint prop", - "MATCH (a:Atom {uid:'A0'})-[e]->{1,8}(b:Atom) WHERE b.element = 'O' RETURN count(b)") - - print("\n=== 2. NULL handling (COD has many missing fields) ===") - # temperature exists on both Structures; add a node missing it to test. - try: - client.new_change() - client.query( - "CREATE (:Structure {cod_id:1000003, r_factor:0.09, formula:'C1 H4', " - "hm_symbol:'P 1'})" - ) - client.query("CHANGE SUBMIT") - client.checkout() - print(" (added Structure 1000003 with NO temperature property)") - except Exception as exc: # noqa: BLE001 - print(f" could not add test node: {' '.join(str(exc).split())[:150]}") - - show(client, "IS NOT NULL on a partially-present prop", - "MATCH (s:Structure) WHERE s.temperature IS NOT NULL RETURN s.cod_id") - show(client, "IS NULL on a partially-present prop", - "MATCH (s:Structure) WHERE s.temperature IS NULL RETURN s.cod_id") - show(client, "plain projection of a missing prop", - "MATCH (s:Structure) RETURN s.cod_id, s.temperature") - show(client, "comparison against a missing prop", - "MATCH (s:Structure) WHERE s.temperature < 200.0 RETURN s.cod_id") - show(client, "sentinel approach: has_temperature flag", - "MATCH (s:Structure) WHERE s.r_factor > 0.0 RETURN count(s)") - - print("\n=== 3. Property indexes ===") - try: - client.new_change() - client.query("CREATE INDEX uid_index FOR (n) ON n.uid") - client.query("CHANGE SUBMIT") - client.checkout() - print(" created index uid_index") - except Exception as exc: # noqa: BLE001 - print(f" FAIL creating index: {' '.join(str(exc).split())[:200]}") - show(client, "db.showIndexes()", "CALL db.showIndexes()") - show(client, "seeded lookup via indexed prop", - "MATCH (a:Atom {uid:'A5'}) RETURN a.element") - - print("\n=== 4. LOAD JSONL (documented APOC shape) ===") - data_dir = os.path.join(TURING_DIR, "data") - os.makedirs(data_dir, exist_ok=True) - path = os.path.join(data_dir, "probe_jsonl.jsonl") - lines = [ - {"type": "node", "id": "0", "labels": ["Atom"], - "properties": {"uid": "J0", "element": "O", "occupancy": 1.0}}, - {"type": "node", "id": "1", "labels": ["Atom"], - "properties": {"uid": "J1", "element": "N", "occupancy": 1.0}}, - {"type": "node", "id": "2", "labels": ["Structure"], - "properties": {"cod_id": 2000001, "formula": "C2 H6"}}, - {"type": "relationship", "id": "0", "label": "CONTACT", - "start": {"id": "0"}, "end": {"id": "1"}, - "properties": {"kind": "hbond", "length": 2.11, "symop": "x,y,z"}}, - {"type": "relationship", "id": "1", "label": "HAS_SITE", - "start": {"id": "2"}, "end": {"id": "0"}, "properties": {}}, - ] - with open(path, "w") as fh: - for obj in lines: - fh.write(json.dumps(obj) + "\n") - print(f" wrote {path}") - - for stmt in [ - "LOAD JSONL 'probe_jsonl.jsonl' AS probe_jsonl_graph", - ]: - try: - client.query(stmt) - print(f" PASS {stmt}") - except Exception as exc: # noqa: BLE001 - print(f" FAIL {stmt}\n {' '.join(str(exc).split())[:250]}") - - try: - client.set_graph("probe_jsonl_graph") - show(client, "loaded graph node count", "MATCH (n) RETURN count(n)") - show(client, "loaded graph edge count", "MATCH ()-[r]->() RETURN count(r)") - show(client, "loaded props survived", "MATCH (a:Atom) RETURN a.uid, a.element") - show(client, "loaded edge props survived", - "MATCH ()-[r:CONTACT]->() RETURN r.kind, r.length, r.symop") - show(client, "labels round-tripped", "CALL db.labels()") - show(client, "edge types round-tripped", "CALL db.edgeTypes()") - except Exception as exc: # noqa: BLE001 - print(f" could not query loaded graph: {' '.join(str(exc).split())[:200]}") - - return 0 - - -if __name__ == "__main__": - sys.exit(main()) diff --git a/tests/probe_gemmi.py b/tests/probe_gemmi.py deleted file mode 100644 index f3bd168..0000000 --- a/tests/probe_gemmi.py +++ /dev/null @@ -1,105 +0,0 @@ -"""Inspect gemmi's small-structure + symmetry API against a real COD CIF. - -Symmetry expansion is where this pipeline is most likely to be silently wrong, -so before writing any of it we check what gemmi actually offers on the installed -version rather than assuming an API shape. -""" - -from __future__ import annotations - -import glob -import sys - -import gemmi - - -def members(obj: object, needle: str = "") -> list[str]: - return sorted( - m for m in dir(obj) - if not m.startswith("_") and (needle.lower() in m.lower() if needle else True) - ) - - -def main() -> int: - print(f"gemmi {gemmi.__version__}\n") - - paths = sorted(glob.glob("data/cif/**/*.cif", recursive=True)) - if not paths: - print("no CIFs found under data/cif -- run ingest.download first") - return 1 - print(f"{len(paths)} CIFs available; using {paths[0]}\n") - - print("=== module-level readers ===") - print(members(gemmi, "small")) - print(members(gemmi, "read")) - - st = gemmi.read_small_structure(paths[0]) - print(f"\n=== SmallStructure ({type(st).__name__}) ===") - print(members(st)) - - print(f"\n name = {st.name}") - print(f" spacegroup_hm = {st.spacegroup_hm!r}") - print(f" cell = {st.cell}") - print(f" n sites = {len(st.sites)}") - - print(f" spacegroup_hall = {st.spacegroup_hall!r}") - print(f" spacegroup_number = {st.spacegroup_number}") - print(f" symops ({len(st.symops)}): {[str(o) for o in st.symops][:4]}") - - sg = st.spacegroup - print(f"\n=== .spacegroup -> {sg} ===") - if sg is not None: - print(members(sg)) - print(f" number={sg.number} hm={sg.hm!r} hall={sg.hall!r} " - f"xhm={sg.xhm()!r} crystal_system={sg.crystal_system_str()!r}") - ops = sg.operations() - oplist = list(ops) - print(f" n operations = {len(oplist)}") - for i, op in enumerate(oplist[:4]): - print(f" op[{i}] = {op.triplet()}") - print(f"\n check_spacegroup: {st.check_spacegroup}") - print(f" determine_and_set_spacegroup: {st.determine_and_set_spacegroup}") - - site = st.sites[0] - print(f"\n=== SmallStructure.Site ===") - print(members(site)) - print(f" label={site.label!r} type_symbol={site.type_symbol!r} " - f"element={site.element} occ={site.occ} u_iso={site.u_iso}") - print(f" fract={site.fract} orth={site.orth(st.cell)}") - - print("\n=== unit-cell expansion helpers ===") - for name in ("get_all_unit_cell_sites", "setup_cell_images", "cell"): - print(f" has {name}: {hasattr(st, name)}") - if hasattr(st, "get_all_unit_cell_sites"): - allsites = st.get_all_unit_cell_sites() - print(f" get_all_unit_cell_sites() -> {len(allsites)} sites " - f"(asymmetric unit has {len(st.sites)})") - print(f" site type: {type(allsites[0]).__name__}, members: {members(allsites[0])}") - - print("\n=== NeighborSearch / ContactSearch ===") - print(" NeighborSearch:", members(gemmi.NeighborSearch)) - print(" ContactSearch :", members(gemmi.ContactSearch)) - - print("\n=== UnitCell ===") - print(members(st.cell, "find")) - print(members(st.cell, "image")) - print(members(st.cell, "orth")) - print(members(st.cell, "fract")) - print(f" volume = {st.cell.volume}") - - print("\n=== Element ===") - el = gemmi.Element("O") - print(members(el)) - print(f" O: atomic_number={el.atomic_number} vdw_r={el.vdw_r} " - f"covalent_r={el.covalent_r} weight={el.weight} is_metal={el.is_metal}") - - print("\n=== small structure -> Structure conversion (for NeighborSearch) ===") - for name in ("make_structure_from_small", "small_to_structure"): - print(f" gemmi.{name}: {hasattr(gemmi, name)}") - print(members(gemmi, "structure")) - - return 0 - - -if __name__ == "__main__": - sys.exit(main()) diff --git a/tests/probe_neighbor.py b/tests/probe_neighbor.py deleted file mode 100644 index 0dd0950..0000000 --- a/tests/probe_neighbor.py +++ /dev/null @@ -1,142 +0,0 @@ -"""Verify periodic neighbour search and, crucially, symmetry-operation recovery. - -The spec is explicit that `symop` -- the symmetry operation plus lattice -translation that generated the neighbouring atom's position -- is the technical -heart of the ingest, because it is what lets a query traverse the infinite -periodic contact network without materialising a supercell. - -So before writing the contact detector we establish, against a real COD entry: - * how to build a NeighborSearch over a SmallStructure - * what find_site_neighbors returns - * how to turn a Mark back into (symmetry op index, lattice translation) - * that the recovered symop actually reproduces the neighbour's position -""" - -from __future__ import annotations - -import glob -import inspect -import sys - -import gemmi - - -def sig(fn: object, name: str) -> None: - try: - print(f" {name}{inspect.signature(fn)}") # type: ignore[arg-type] - except (TypeError, ValueError): - doc = (getattr(fn, "__doc__", "") or "").strip().splitlines() - print(f" {name}: {doc[0] if doc else ''}") - - -def main() -> int: - paths = sorted(glob.glob("data/cif/**/*.cif", recursive=True)) - if not paths: - print("no CIFs; run ingest.download first") - return 1 - - # pick a structure with a decent number of sites and a non-trivial group - chosen = None - for p in paths[:400]: - try: - s = gemmi.read_small_structure(p) - except Exception: # noqa: BLE001 - continue - if s.spacegroup is not None and 10 <= len(s.sites) <= 60 and s.spacegroup.number > 1: - chosen = (p, s) - break - if chosen is None: - print("no suitable structure found") - return 1 - path, st = chosen - st.setup_cell_images() - print(f"using {path}") - print(f" {st.name} {st.spacegroup_hm} (#{st.spacegroup.number}) " - f"{len(st.sites)} sites cell={st.cell}\n") - - print("=== constructor / method signatures ===") - sig(gemmi.NeighborSearch.__init__, "NeighborSearch.__init__") - for m in ("add_site", "populate", "find_site_neighbors", "find_atoms", - "get_image_transformation", "find_neighbors"): - sig(getattr(gemmi.NeighborSearch, m, None), f"NeighborSearch.{m}") - print() - sig(getattr(gemmi.UnitCell, "find_nearest_pbc_image", None), - "UnitCell.find_nearest_pbc_image") - sig(getattr(gemmi.UnitCell, "find_nearest_pbc_position", None), - "UnitCell.find_nearest_pbc_position") - - print("\n=== build the search ===") - try: - ns = gemmi.NeighborSearch(st, 5.0) - print(" gemmi.NeighborSearch(small_structure, max_radius) OK") - except Exception as exc: # noqa: BLE001 - print(f" small-structure ctor failed: {exc}") - return 1 - ns.populate() - print(f" populated. grid_cell={ns.grid_cell} radius_specified={ns.radius_specified}") - - print("\n=== Mark structure ===") - print(" Mark members:", sorted(m for m in dir(gemmi.NeighborSearch.Mark) - if not m.startswith("_"))) - - print("\n=== find_site_neighbors on site 0 ===") - site0 = st.sites[0] - marks = ns.find_site_neighbors(site0, min_dist=0.1, max_dist=4.0) - print(f" site {site0.label} ({site0.type_symbol}) -> {len(marks)} neighbours") - - all_sites = st.get_all_unit_cell_sites() - print(f" unit cell has {len(all_sites)} sites (asym {len(st.sites)})") - - ops = list(st.spacegroup.operations()) - print(f" spacegroup has {len(ops)} operations") - - print("\n=== recovering the symmetry operation per contact ===") - pos0 = st.cell.orthogonalize(site0.fract) - shown = 0 - for mark in marks: - if shown >= 6: - break - # image_idx 0 is the identity; >0 indexes cell images built by - # setup_cell_images(). to_site maps the mark back to a unit-cell site. - img = getattr(mark, "image_idx", None) - try: - ref = mark.to_site(all_sites) - ref_label = ref.label - ref_el = ref.type_symbol - except Exception as exc: # noqa: BLE001 - ref_label = f"" - ref_el = "?" - # the actual position of this periodic image - mpos = mark.pos - d = pos0.dist(mpos) - # what transformation took the reference site to this image? - try: - ft = ns.get_image_transformation(img) if img is not None else None - tr = f"mat={[round(v,3) for row in ft.mat.tolist() for v in row]} vec={[round(v,4) for v in ft.vec.tolist()]}" if ft else "identity" - except Exception as exc: # noqa: BLE001 - tr = f"" - print(f" {site0.label}...{ref_label:<6} ({ref_el:<2}) d={d:5.3f} " - f"image_idx={img}\n {tr}") - shown += 1 - - print("\n=== ContactSearch (the higher-level API) ===") - for m in ("search_radius", "ignore", "twice", "min_occupancy", - "special_pos_cutoff_sq", "setup_atomic_radii", "set_radius", - "get_radius", "find_contacts"): - sig(getattr(gemmi.ContactSearch, m, None), f"ContactSearch.{m}") - print(" Ignore enum:", sorted(m for m in dir(gemmi.ContactSearch.Ignore) - if not m.startswith("_"))) - print(" Result members:", sorted(m for m in dir(gemmi.ContactSearch.Result) - if not m.startswith("_"))) - - print("\n=== gemmi's built-in radii (to decide whether to ship our own) ===") - for symbol in ("H", "C", "N", "O", "F", "S", "Cl", "Br", "I", "Cu", "Zn", "Fe"): - el = gemmi.Element(symbol) - print(f" {symbol:<3} Z={el.atomic_number:<3} covalent_r={el.covalent_r:.3f} " - f"vdw_r={el.vdw_r:.3f} is_metal={el.is_metal}") - - return 0 - - -if __name__ == "__main__": - sys.exit(main()) diff --git a/tests/probe_turingdb_syntax.py b/tests/probe_turingdb_syntax.py deleted file mode 100644 index af55a83..0000000 --- a/tests/probe_turingdb_syntax.py +++ /dev/null @@ -1,194 +0,0 @@ -"""Probe TuringDB's OWN documented syntax (not Neo4j's) for the demo's needs. - -probe_critical.py deliberately tests the Neo4j spellings from the project spec. -This file tests the forms docs.turingdb.ai actually documents, so we learn what -the demo CAN express rather than only what it can't: - - * variable length is a POSTFIX QUANTIFIER -[e]->{2,4}(b) not -[e*2..4]->(b) - * shortestPath is a bare STATEMENT with 5 args, not a pattern function - * aggregates are count()/avg() only, and only as a SINGLE return item - * property indexes exist and matter for seeded lookups - -Usage: - TURING_HOST=http://localhost:6691 ./.venv/bin/python tests/probe_turingdb_syntax.py -""" - -from __future__ import annotations - -import os -import sys - -from turingdb import TuringDB, TuringDBException - -HOST = os.environ.get("TURING_HOST", "http://localhost:6691") -GRAPH = os.environ.get("PROBE_GRAPH", "probe_critical") # reuse the built fixture - -# (group, label, cypher, why) -TESTS: list[tuple[str, str, str, str]] = [ - # ---------------------------------------------------------------- Q2 --- - # The headline claim. Docs give postfix +, * and {m,n}. - ("quantifier", "bare edge, one-or-more -[e]->+", - "MATCH (a:Atom)-[e]->+(b:Atom) RETURN count(b)", "Q2"), - ("quantifier", "bare edge, zero-or-more -[e]->*", - "MATCH (a:Atom)-[e]->*(b:Atom) RETURN count(b)", "Q2"), - ("quantifier", "bare edge, range -[e]->{2,4}", - "MATCH (a:Atom)-[e]->{2,4}(b:Atom) RETURN count(b)", "Q2"), - ("quantifier", "TYPED edge + range -[e:CONTACT]->{1,8}", - "MATCH (a:Atom)-[e:CONTACT]->{1,8}(b:Atom) RETURN count(b)", - "Q2 restricted to contacts - docs do NOT confirm this combination"), - ("quantifier", "typed anon edge + range -[:CONTACT]->{1,8}", - "MATCH (a:Atom)-[:CONTACT]->{1,8}(b:Atom) RETURN count(b)", "Q2"), - ("quantifier", "SEEDED typed + range (the real Q2)", - "MATCH (a:Atom {uid:'A0'})-[e:CONTACT]->{1,8}(b:Atom) RETURN count(b)", - "Q2 exactly as the demo needs it"), - ("quantifier", "exact depth {8,8}", - "MATCH (a:Atom {uid:'A0'})-[e:CONTACT]->{8,8}(b:Atom) RETURN count(b)", - "Q2 growth curve, one query per depth"), - ("quantifier", "exact depth {3}", - "MATCH (a:Atom {uid:'A0'})-[e:CONTACT]->{3}(b:Atom) RETURN count(b)", - "Q2 shorthand"), - ("quantifier", "quantifier + WHERE on rel prop", - "MATCH (a:Atom {uid:'A0'})-[e:CONTACT]->{1,8}(b:Atom) WHERE e.kind = 'hbond' " - "RETURN count(b)", "Q2 hbond-only network"), - ("quantifier", "quantifier + INLINE rel prop", - "MATCH (a:Atom {uid:'A0'})-[e:CONTACT {kind:'hbond'}]->{1,8}(b:Atom) RETURN count(b)", - "Q2 hbond-only network"), - ("quantifier", "undirected quantifier -[e]-{1,3}", - "MATCH (a:Atom {uid:'A0'})-[e:CONTACT]-{1,3}(b:Atom) RETURN count(b)", - "packing network ignores donor/acceptor direction"), - ("quantifier", "reverse quantifier <-[e]-{1,3}", - "MATCH (a:Atom {uid:'A0'})<-[e:CONTACT]-{1,3}(b:Atom) RETURN count(b)", "Q2"), - ("quantifier", "quantifier reaching 12 hops", - "MATCH (a:Atom {uid:'A0'})-[e:CONTACT]->{1,12}(b:Atom) RETURN count(b)", - "depth headroom beyond 8"), - - # ------------------------------------------------------------ shortest --- - ("shortestpath", "shortestPath statement form", - "MATCH (n:Atom {uid:'A0'}), (m:Atom {uid:'A5'}) " - "shortestPath(n, m, length, dist, path) RETURN dist, path", - "path between two sites"), - - # ---------------------------------------------------------- aggregates --- - ("aggregate", "count(*) single item", - "MATCH (a:Atom) RETURN count(*)", "baseline"), - ("aggregate", "count + avg combined in ONE item", - "MATCH ()-[e:CONTACT]->() RETURN count(e) + avg(e.length)", - "docs say aggregates may be combined within a single return item"), - ("aggregate", "avg alone", "MATCH ()-[e:CONTACT]->() RETURN avg(e.length)", "reporting"), - ("aggregate", "count with WHERE = per-type frequency (Q4 workaround)", - "MATCH (f:Fragment)<-[:HAS_FRAGMENT]-(c:Component) " - "WHERE f.fragment_type = 'carboxylic_acid' RETURN count(c)", - "Q4 motif table: one query per fragment type instead of GROUP BY"), - ("aggregate", "count over a 4-hop join (Q4 shape)", - "MATCH (fa:Fragment)<-[:HAS_FRAGMENT]-(c:Component)<-[:CONTAINS_COMPONENT]-(s:Structure) " - "WHERE fa.fragment_type = 'carboxylic_acid' RETURN count(s)", - "Q4 per-motif count"), - - # --------------------------------------------------------------- joins --- - ("join", "comma-joined 3 patterns sharing vars", - "MATCH (f:Fragment)<-[:HAS_FRAGMENT]-(c:Component), " - "(c)<-[:CONTAINS_COMPONENT]-(s:Structure), (s)-[:HAS_SITE]->(a:Atom) " - "RETURN s.cod_id, f.fragment_type, a.uid", "Q1/Q3 need multi-pattern joins"), - ("join", "self-join on shared node with <>", - "MATCH (c1:Component)-[:HAS_FRAGMENT]->(f:Fragment)<-[:HAS_FRAGMENT]-(c2:Component) " - "WHERE c1.inchikey <> c2.inchikey RETURN c1.name, c2.name", - "Q3 chemical-similarity step"), - ("join", "Q3 full 5-hop chain, no aggregate", - "MATCH (api:Component {inchikey:'KDY'})-[:HAS_FRAGMENT]->(f:Fragment)" - "<-[:HAS_FRAGMENT]-(sim:Component), " - "(sim)<-[:CONTAINS_COMPONENT]-(st:Structure), " - "(st)-[:CONTAINS_COMPONENT]->(cand:Component) " - "RETURN cand.inchikey, cand.name", "Q3 core traversal"), - ("join", "Q1 cycle: two rels between the same pair", - "MATCH (a1:Atom)-[h1:CONTACT]->(a2:Atom), (a2)-[h2:CONTACT]->(a1) " - "RETURN a1.uid, a2.uid", "Q1 acid dimer - reciprocal pair"), - ("join", "Q1 cycle via distinct endpoint vars", - "MATCH (a1:Atom)-[h1:CONTACT]->(a2:Atom), (a3:Atom)-[h2:CONTACT]->(a4:Atom) " - "WHERE a1.uid = a4.uid AND a2.uid = a3.uid RETURN a1.uid, a2.uid", - "Q1 workaround if the direct cycle will not plan"), - - # -------------------------------------------------------------- indexes --- - ("index", "CALL db.showIndexes()", "CALL db.showIndexes()", "seed-lookup speed"), - ("index", "CALL db.procedures()", "CALL db.procedures()", "discover real procedure list"), - - # ------------------------------------------------------- introspection --- - ("introspection", "CALL db.labels()", "CALL db.labels()", "schema check"), - ("introspection", "CALL db.edgeTypes()", "CALL db.edgeTypes()", "schema check"), - ("introspection", "CALL db.propertyTypes()", "CALL db.propertyTypes()", "schema check"), - ("introspection", "CALL db.history()", "CALL db.history()", "Stage 4 commit log"), - ("introspection", "CALL ... YIELD ... RETURN *", - "CALL db.propertyTypes() YIELD propertyType, valueType RETURN *", "docs example"), - ("introspection", "labels(n) + edgeType(e)", - "MATCH (n)-[e]->(m) RETURN labels(n), edgeType(e), labels(m) LIMIT 5", - "UI node/edge typing"), - ("introspection", "WHERE n = ", - "MATCH (n) WHERE n = 1 RETURN labels(n)", "seed by internal id"), - ("introspection", "CHANGE LIST", "CHANGE LIST", "Stage 4"), - - # --------------------------------------------------------------- misc ---- - ("misc", "RETURN with AS alias", - "MATCH (s:Structure) RETURN s.cod_id AS cod_id", "API payload shaping"), - ("misc", "ORDER BY projected property", - "MATCH ()-[e:CONTACT]->() RETURN e.length ORDER BY e.length LIMIT 5", "Q1 ordering"), - ("misc", "arithmetic in RETURN", - "MATCH (s:Structure) RETURN s.r_factor * 100.0", "reporting"), - ("misc", "IS NULL", - "MATCH (s:Structure) WHERE s.pressure IS NULL RETURN s.cod_id", - "COD missing-field handling"), - ("misc", "IS NOT NULL", - "MATCH (s:Structure) WHERE s.temperature IS NOT NULL RETURN s.cod_id", - "COD missing-field handling"), - ("misc", "list literal RETURN", - "RETURN [1,2,3] AS nums", "docs say literal elements only"), -] - - -def main() -> int: - print(f"connecting to {HOST}, graph={GRAPH}") - client = TuringDB(host=HOST) - try: - client.query(f"LOAD GRAPH {GRAPH}") - except Exception: # noqa: BLE001 - pass - client.set_graph(GRAPH) - n = client.query("MATCH (n) RETURN count(n)").iloc[0, 0] - print(f"fixture has {n} nodes\n") - - width = max(len(t[1]) for t in TESTS) - passes: list[tuple[str, str, str]] = [] - fails: list[tuple[str, str, str, str]] = [] - current = "" - for group, label, cypher, why in TESTS: - if group != current: - print(f"\n--- {group} ---") - current = group - try: - df = client.query(cypher) - except Exception as exc: # noqa: BLE001 - msg = " ".join(str(exc).split()) - print(f" FAIL {label:<{width}} {msg[:110]}") - fails.append((group, label, why, msg[:300])) - else: - val = "" - try: - if len(df) and len(df.columns) == 1: - val = f" = {df.iloc[0, 0]}" - elif len(df): - val = f" first={list(df.iloc[0])[:4]}" - except Exception: # noqa: BLE001 - pass - print(f" PASS {label:<{width}} rows={len(df)}{val}") - passes.append((group, label, why)) - - print(f"\n{'=' * 70}\n{len(passes)}/{len(TESTS)} supported") - if fails: - print("\nUNSUPPORTED:") - for group, label, why, msg in fails: - print(f" [{group}] {label}") - print(f" needed for: {why}") - print(f" error : {msg[:180]}") - return 0 - - -if __name__ == "__main__": - sys.exit(main()) diff --git a/tests/validate_chemistry.py b/tests/validate_chemistry.py index 2babd72..04d8026 100644 --- a/tests/validate_chemistry.py +++ b/tests/validate_chemistry.py @@ -18,6 +18,10 @@ import gemmi +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) + from ingest import chemistry as ch from ingest.structure import build_atoms, detect_bonds, find_components, read_structure @@ -60,6 +64,10 @@ def _is_metal(symbol: str) -> bool: def main() -> int: n = int(sys.argv[1]) if len(sys.argv) > 1 else 300 paths = sorted(glob.glob("data/cif/**/*.cif", recursive=True)) + if not paths: + print("no CIFs under data/cif -- run `python -m ingest.download " + "--dataset slice` first", file=sys.stderr) + return 1 sample = random.Random(11).sample(paths, min(n, len(paths))) valence: dict[str, Counter[int]] = defaultdict(Counter) @@ -121,6 +129,14 @@ def main() -> int: f"({100 * isolated_heavy / max(1, total_heavy):.2f}%) " f"[excludes halide counter-ions]") print("\nvalence in PURELY ORGANIC environments (no metal neighbour):") + # A validator that analysed nothing must not report PASS. Without this, + # every element is skipped on an empty corpus, `ok` is never cleared, and + # the script exits 0 having checked zero structures -- which reads exactly + # like a clean run. + if not any(valence.get(el) for el in EXPECTED): + print(" no atoms analysed -- run `python -m ingest.download --dataset slice` " + "first", file=sys.stderr) + return 1 ok = True for el in ("H", "C", "N", "O", "F", "S", "Cl", "Br", "I"): counts = valence.get(el) diff --git a/tests/validate_components.py b/tests/validate_components.py index d349f79..e7bafa1 100644 --- a/tests/validate_components.py +++ b/tests/validate_components.py @@ -14,7 +14,10 @@ import gemmi -from ingest import chemistry as ch +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) + from ingest.structure import ( build_atoms, component_instances, @@ -44,6 +47,10 @@ def main() -> int: check_op_api() paths = sorted(glob.glob("data/cif/**/*.cif", recursive=True)) + if not paths: + print("no CIFs under data/cif -- run `python -m ingest.download " + "--dataset slice` first", file=sys.stderr) + return 1 sample = random.Random(7).sample(paths, min(n * 6, len(paths))) shown = 0 diff --git a/tests/validate_symmetry.py b/tests/validate_symmetry.py index 825dd5e..222d7ae 100644 --- a/tests/validate_symmetry.py +++ b/tests/validate_symmetry.py @@ -18,6 +18,10 @@ import gemmi +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) + from ingest import chemistry as ch from ingest.symmetry import resolve_symop, spacegroup_of, verify_symop diff --git a/visualizer/src/components/viewer/crystal/chat.tsx b/visualizer/src/components/viewer/crystal/chat.tsx index b2a041e..5ef0cd9 100644 --- a/visualizer/src/components/viewer/crystal/chat.tsx +++ b/visualizer/src/components/viewer/crystal/chat.tsx @@ -11,16 +11,20 @@ import { type FC, useCallback, useEffect, useRef, useState } from 'react' import { Icon } from '@blueprintjs/core' import { useAppStore, useVisStore } from '@/stores' +import { isCypherQuery } from '@/utils/is-cypher' import { useTuringContext } from '@turingcanvas' import { CRYSTAL_GRAPHS, INVESTIGATIONS, + SHELVES, type Investigation, type Snapshot, VERSIONED_GRAPH, loadLedger, loadSubgraph, routeQuestion, + runRawCypher, + type RawResult, } from './model' import { paintSubgraph } from './canvas' @@ -38,6 +42,8 @@ interface Turn { shown: number // words revealed so far, across the whole answer /** set when the investigation can be replayed across commits */ inv?: Investigation + /** set when the user typed Cypher instead of picking an investigation */ + raw?: RawResult } const GREETING = 'Good afternoon, what would you like to investigate' @@ -136,6 +142,55 @@ export const CrystalChat: FC = () => { [graphName, turing] ) + /** + * Execute Cypher the user typed. Goes through this studio's own painter + * rather than the stock toolbar's pipeline, which is disabled on these + * graphs because it would fetch all 4.7M atoms and overwrite the canvas. + */ + const runRaw = useCallback( + async (query: string) => { + if (!graphName) return + if (timerRef.current) window.clearInterval(timerRef.current) + setErr('') + setShowCypher(true) + setPhase('querying') + setMetrics([]) + setActiveCommit(null) + setTurn({ + question: 'Your query', + cypher: query, + nodes: 0, + edges: 0, + ms: 0, + paragraphs: [], + shown: 0, + }) + try { + const res = await runRawCypher(graphName, query) + if (res.kind === 'graph' && res.sub) { + setPhase('painting') + paintSubgraph(turing, res.sub) + } + setTurn((t) => + t + ? { + ...t, + nodes: res.sub?.nodes.length ?? 0, + edges: res.sub?.edges.length ?? 0, + ms: res.ms, + raw: res, + } + : t + ) + setPhase('done') + } catch (e) { + setErr(e instanceof Error ? e.message : String(e)) + setPhase('idle') + } + }, + [graphName, turing] + ) + /** Re-run the current investigation against one commit and repaint. */ const replayAt = useCallback( async (snap: Snapshot | null) => { @@ -156,13 +211,27 @@ export const CrystalChat: FC = () => { t ? { ...t, nodes: sub.nodes.length, edges: sub.edges.length } : t ) setMetrics(await inv.timeline.metrics(snap ? snap.commit : '')) + + // The prose has to move with the chips, or the paragraph describes + // HEAD while the chips describe the selected commit. + if (snap && inv.timeline.answerAt) { + const paras = await inv.timeline.answerAt(snap.commit, snap.tag) + const words = paras.reduce((n, x) => n + x.split(' ').length, 0) + setTurn((t) => (t ? { ...t, paragraphs: paras, shown: words } : t)) + } else if (!snap) { + // back to HEAD: restore the investigation's own narrative + const target = inv.graph ?? graphName ?? '' + const paras = await inv.answer(target, sub) + const words = paras.reduce((n, x) => n + x.split(' ').length, 0) + setTurn((t) => (t ? { ...t, paragraphs: paras, shown: words } : t)) + } } catch (e) { setErr(e instanceof Error ? e.message : String(e)) } finally { setTlBusy(false) } }, - [turn?.inv, turn?.cypher, turing] + [turn?.inv, turn?.cypher, turing, graphName] ) useEffect(() => { @@ -175,7 +244,11 @@ export const CrystalChat: FC = () => { const q = input.trim() if (!q) return setInput('') - run(routeQuestion(q), q) + // Anything starting MATCH / CALL / RETURN / ... is meant as Cypher; the + // rest is routed to the nearest investigation by word overlap. Same + // detection the stock toolbar uses, so the two behave consistently. + if (isCypherQuery(q)) runRaw(q) + else run(routeQuestion(q), q) } // reveal the answer progressively across paragraphs @@ -285,31 +358,78 @@ export const CrystalChat: FC = () => { marginBottom: 9, }} > - Suggested - -
- {INVESTIGATIONS.map((inv) => ( - - ))} + Query library
+ {/* Grouped rather than one flat list: a CCDC room splits into + crystallographers, solid-form/design people and the database + team, and each shelf is aimed at one of them. Every entry + carries a one-line statement of what it DEMONSTRATES, since + the prompt alone does not say why the question is hard. */} + {SHELVES.map((shelf) => { + const items = INVESTIGATIONS.filter((i) => i.shelf === shelf.id) + if (!items.length) return null + return ( +
+
+ {shelf.title} +
+
+ {shelf.note} +
+
+ {items.map((inv) => ( + + ))} +
+
+ ) + })} )} @@ -386,6 +506,115 @@ export const CrystalChat: FC = () => { )} + {/* Put the query in the box so it can be edited and re-run. + Reading the Cypher is worth something; changing one + predicate and watching the canvas move is worth much more + to an audience that writes queries for a living. */} + {showCypher && ( + + )} + + {turn.raw?.kind === 'table' && ( +
+ + + + {turn.raw.columns?.map((c) => ( + + ))} + + + + {turn.raw.rows?.map((r, i) => ( + + {r.map((v, j) => ( + + ))} + + ))} + +
+ {c} +
+ {v} +
+ {turn.raw.truncated && ( +
+ first 40 rows shown +
+ )} +
+ )} + + {turn.raw?.kind === 'graph' && turn.nodes === 0 && ( +
+ The query ran and returned no rows. +
+ )} + {busy && (
{ onKeyDown={(e) => { if (e.key === 'Enter') submit() }} - placeholder="Ask about the corpus…" + placeholder="Ask a question, or type Cypher…" style={{ flex: 1, background: '#12151a', diff --git a/visualizer/src/components/viewer/crystal/index.tsx b/visualizer/src/components/viewer/crystal/index.tsx index 7921df4..c66fe9c 100644 --- a/visualizer/src/components/viewer/crystal/index.tsx +++ b/visualizer/src/components/viewer/crystal/index.tsx @@ -223,8 +223,18 @@ export const CrystalPanel: FC = () => { lineHeight: 1.6, }} > + {/* Three populations, shown as three, because a blended contact + count is the number that gets challenged first. */}
- {n(overview.hbonds)} hydrogen bonds,{' '} + + {n(overview.hbonds - overview.hbondsInferred)} + {' '} + strong H-bonds (H located),{' '} + {n(overview.hbondsInferred)}{' '} + inferred +
+
+ {n(overview.weak)} weak C–H···A,{' '} {n(overview.halogen)} halogen bonds
@@ -235,12 +245,19 @@ export const CrystalPanel: FC = () => { at{' '} {overview.meanAngle.toFixed(1)}° + {' '} + + (strong, located-H population only)
- {n(overview.hbondsInferred)} hydrogen bonds are h_inferred — no refined H - position, so a heavy-atom cutoff was used. Kept flagged, never mixed in - silently. + {overview.hbonds + ? ((100 * overview.hbondsInferred) / overview.hbonds).toFixed(0) + : '0'} + % of the hydrogen-bond layer is h_inferred — no refined H position, so a + heavy-atom D···A cutoff was used. Only N, O and S are admitted as inferred + donors; {n(overview.closeContact)} contacts with an implausible donor are + stored as close_contact and excluded from every H-bond statistic.
)} diff --git a/visualizer/src/components/viewer/crystal/model.ts b/visualizer/src/components/viewer/crystal/model.ts index 82941c7..9ad0da6 100644 --- a/visualizer/src/components/viewer/crystal/model.ts +++ b/visualizer/src/components/viewer/crystal/model.ts @@ -191,6 +191,121 @@ export async function loadSubgraph( return { nodes: [...nodes.values()], edges: [...edges.values()] } } +/** + * Run a Cypher string the user typed, rather than one of the canned + * investigations. + * + * The stock visualizer toolbar is deliberately hidden on the crystallography + * graphs (see pages/viewer.tsx): its pipeline fetches EVERY node and edge in + * the graph, which is hopeless at 4.7M atoms, and its next sync would wipe the + * subgraph this studio paints by hand. That is a good reason to hide THAT + * toolbar and a bad reason to have no Cypher box at all -- being able to read + * and edit the query is most of what makes the demo credible to an audience + * that writes queries for a living. + * + * So raw Cypher goes through the studio's own painter instead. Two shapes come + * back: + * + * - a query projecting the 8-column subgraph shape + * (`a, labels(a), , e, edgeType(e), b, labels(b), `) paints the + * canvas, exactly as an investigation does; + * - anything else -- a count, an average, a flat projection -- comes back as + * columns to tabulate, because refusing to show a scalar would make the box + * useless for precisely the checks someone wants to run live. + */ +export interface RawResult { + kind: 'graph' | 'table' + sub?: Subgraph + columns?: string[] + rows?: string[][] + truncated?: boolean + ms: number +} + +/** Rows to show before truncating a tabular result. */ +const RAW_ROW_LIMIT = 40 + +export async function runRawCypher( + graph: string, + query: string +): Promise { + const t0 = performance.now() + const res = await fetch(`/api/query?graph=${encodeURIComponent(graph)}`, { + method: 'POST', + headers: { 'Content-Type': 'text/plain' }, + body: query, + }) + const json = await res.json() + // Failures arrive as HTTP 200 with a non-null top-level `error`. Surface the + // engine's own wording -- PARSE_ERROR vs ANALYZE_ERROR vs PLAN_ERROR means + // something quite specific here, and paraphrasing it would hide that. + if (json.error) { + throw new Error( + `${json.error}${json.error_details ? `: ${json.error_details}` : ''}` + ) + } + const ms = performance.now() - t0 + const names: string[] = json.header?.column_names ?? [] + const cols = columns(json.data as unknown[][]) + + const looksLikeSubgraph = + names.length === 8 && + /^labels\(/i.test(names[1] ?? '') && + /^(edgeType|type)\(/i.test(names[4] ?? '') && + /^labels\(/i.test(names[6] ?? '') + + if (looksLikeSubgraph) { + const nodes = new Map() + const edges = new Map() + const n = cols[0]?.length ?? 0 + for (let i = 0; i < n; i++) { + const aId = num(cols[0]?.[i]) + const bId = num(cols[5]?.[i]) + const eId = num(cols[3]?.[i]) + if (!nodes.has(aId)) + nodes.set(aId, { + id: aId, + label: str(cols[1]?.[i]), + name: str(cols[2]?.[i]) || str(cols[1]?.[i]), + }) + if (!nodes.has(bId)) + nodes.set(bId, { + id: bId, + label: str(cols[6]?.[i]), + name: str(cols[7]?.[i]) || str(cols[6]?.[i]), + }) + if (!edges.has(eId)) + edges.set(eId, { id: eId, src: aId, tgt: bId, type: str(cols[4]?.[i]) }) + } + return { + kind: 'graph', + sub: { nodes: [...nodes.values()], edges: [...edges.values()] }, + ms, + } + } + + const rowCount = cols[0]?.length ?? 0 + const shown = Math.min(rowCount, RAW_ROW_LIMIT) + const rows: string[][] = [] + for (let i = 0; i < shown; i++) { + rows.push(names.map((_, c) => { + const v = cols[c]?.[i] + if (v === null || v === undefined) return '--' + if (typeof v === 'number') + return Number.isInteger(v) ? v.toLocaleString() : v.toFixed(4) + if (typeof v === 'boolean') return v ? 'true' : 'false' + return String(v) + })) + } + return { + kind: 'table', + columns: names, + rows, + truncated: rowCount > shown, + ms, + } +} + // --------------------------------------------------------------------------- // the commit ledger // --------------------------------------------------------------------------- @@ -247,6 +362,10 @@ export interface Overview { polymeric: number meanHA: number meanAngle: number + /** weak C-H...A contacts: a separate population, never merged with hbond */ + weak: number + /** inferred contacts with an implausible donor, excluded from H-bond stats */ + closeContact: number } export async function loadOverview(graph: string): Promise { @@ -254,25 +373,34 @@ export async function loadOverview(graph: string): Promise { const [ structures, atoms, components, contacts, hbonds, hbondsInferred, halogen, bonds, fragments, spaceGroups, withInchikey, polymeric, meanHA, meanAngle, + weak, closeContact, ] = await Promise.all([ q('MATCH (n:Structure) RETURN count(n)'), q('MATCH (n:Atom) RETURN count(n)'), q('MATCH (n:Component) RETURN count(n)'), q('MATCH ()-[r:CONTACT]->() RETURN count(r)'), q("MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' RETURN count(r)"), - q('MATCH ()-[r:CONTACT]->() WHERE r.h_inferred = true RETURN count(r)'), + q("MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' AND r.h_inferred = true RETURN count(r)"), q("MATCH ()-[r:CONTACT]->() WHERE r.kind = 'halogen' RETURN count(r)"), q('MATCH ()-[r:BONDED_TO]->() RETURN count(r)'), q('MATCH (n:Fragment) RETURN count(n)'), q('MATCH (n:SpaceGroup) RETURN count(n)'), q('MATCH (c:Component) WHERE c.has_inchikey = true RETURN count(c)'), q('MATCH (c:Component) WHERE c.is_polymeric = true RETURN count(c)'), - q('MATCH ()-[r:CONTACT]->() WHERE r.h_inferred = false RETURN avg(r.length)'), - q('MATCH ()-[r:CONTACT]->() WHERE r.h_inferred = false RETURN avg(r.angle)'), + // Scoped to kind='hbond' AND a located hydrogen. Filtering on h_inferred + // alone silently averaged H...A together with the halogen bonds' X...A -- + // and now also with the weak C-H...A population, whose mean sits 0.6 A + // further out. The published 2.09 A / 160.5 deg came from the blended + // version; the located-H hydrogen-bond value is 2.048 A / 160.4 deg. + q("MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' AND r.h_inferred = false RETURN avg(r.length)"), + q("MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' AND r.h_inferred = false RETURN avg(r.angle)"), + q("MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond_weak' RETURN count(r)"), + q("MATCH ()-[r:CONTACT]->() WHERE r.kind = 'close_contact' RETURN count(r)"), ]) return { structures, atoms, components, contacts, hbonds, hbondsInferred, halogen, bonds, fragments, spaceGroups, withInchikey, polymeric, meanHA, meanAngle, + weak, closeContact, } } @@ -363,8 +491,38 @@ export async function loadNodeDetail( // the investigations offered in the chat panel // --------------------------------------------------------------------------- +/** + * The three shelves the query library is arranged on. A CCDC audience splits + * cleanly this way: the crystallographers care about the first, the solid-form + * and drug-design people about the second, and the database team about the + * third. + */ +export type Shelf = 'packing' | 'solidform' | 'provenance' + +export const SHELVES: { id: Shelf; title: string; note: string }[] = [ + { + id: 'packing', + title: 'Packing & motifs', + note: 'The contact network itself -- the edge set neither COD nor the CSD stores.', + }, + { + id: 'solidform', + title: 'Solid form & design', + note: 'Salt vs co-crystal, coformers, polymorphs: the questions that need identity to be cross-structure.', + }, + { + id: 'provenance', + title: 'Provenance & versioning', + note: 'Git-like commits over the corpus: reproduce a published statistic, or find what a correction invalidates.', + }, +] + export interface Investigation { id: string + /** Which shelf of the library this sits on. */ + shelf: Shelf + /** One line, shown under the prompt: what this query DEMONSTRATES. */ + blurb: string prompt: string cypher: string /** Run against a different graph than the one the studio is viewing. */ @@ -382,6 +540,14 @@ export interface Investigation { timeline?: { cypher?: string metrics: (commit: string) => Promise<{ label: string; value: string }[]> + /** + * Narrative for ONE commit. + * + * Without it the narrative would keep describing HEAD while the metric + * chips describe the selected year -- two different corpora side by side + * in the same panel. + */ + answerAt?: (commit: string, tag: string) => Promise } /** Builds the narrative answer from live follow-up queries. */ answer: (graph: string, sub: Subgraph) => Promise @@ -394,9 +560,26 @@ const pct = (a: number, b: number) => (b ? ((100 * a) / b).toFixed(1) : '0.0') /** A structure that actually has a rich contact network, for the packing view. */ export const PACKING_SEED_COD = 2229029 +/** + * Prefix of the acid-dimer counting query, up to the `is_involution` test. + * + * Built once because four follow-ups differ only in that boolean and in what + * they aggregate. Written as ONE connected chain starting at Fragment (16 + * nodes) rather than as comma-separated patterns joined in the WHERE -- that + * rewrite alone is worth roughly 20x in this engine, which no index replaces. + */ +const DIMER_COUNT = + "MATCH (f1:Fragment)<-[:IN_FRAGMENT]-(a1:Atom)" + + "-[h:CONTACT]->(a2:Atom)-[:IN_FRAGMENT]->(f2:Fragment) " + + "WHERE f1.fragment_type = 'carboxylic_acid' AND f2.fragment_type = 'carboxylic_acid' " + + "AND h.kind = 'hbond' AND h.h_inferred = false AND h.is_involution = " + export const INVESTIGATIONS: Investigation[] = [ { id: 'packing', + shelf: 'packing', + blurb: + 'Every contact edge carries the symmetry operation that generated the neighbour, so the periodic network is traversable without building a supercell.', prompt: `Show the hydrogen-bond network of one structure (COD ${PACKING_SEED_COD})`, cypher: `MATCH (a:Atom)-[e:CONTACT]->(b:Atom) WHERE a.cod_id = ${PACKING_SEED_COD} @@ -406,7 +589,7 @@ RETURN a, labels(a), a.label, e, edgeType(e), b, labels(b), b.label`, const cod = PACKING_SEED_COD const [nH, nInf, nHal, mLen, mAng, nComp, nAtoms] = await Promise.all([ scalar(g, `MATCH (a:Atom)-[r:CONTACT]->(b:Atom) WHERE a.cod_id = ${cod} AND r.kind = 'hbond' RETURN count(r)`), - scalar(g, `MATCH (a:Atom)-[r:CONTACT]->(b:Atom) WHERE a.cod_id = ${cod} AND r.h_inferred = true RETURN count(r)`), + scalar(g, `MATCH (a:Atom)-[r:CONTACT]->(b:Atom) WHERE a.cod_id = ${cod} AND r.kind = 'hbond' AND r.h_inferred = true RETURN count(r)`), scalar(g, `MATCH (a:Atom)-[r:CONTACT]->(b:Atom) WHERE a.cod_id = ${cod} AND r.kind = 'halogen' RETURN count(r)`), scalar(g, `MATCH (a:Atom)-[r:CONTACT]->(b:Atom) WHERE a.cod_id = ${cod} AND r.h_inferred = false RETURN avg(r.length)`), scalar(g, `MATCH (a:Atom)-[r:CONTACT]->(b:Atom) WHERE a.cod_id = ${cod} AND r.h_inferred = false RETURN avg(r.angle)`), @@ -423,6 +606,9 @@ RETURN a, labels(a), a.label, e, edgeType(e), b, labels(b), b.label`, }, { id: 'recurrence', + shelf: 'solidform', + blurb: + 'One molecule appearing in forty structures is one node with forty edges, not forty rows joined on a string.', prompt: 'Which molecules recur across different structures?', cypher: `MATCH (s:Structure)-[e:CONTAINS_COMPONENT]->(c:Component) RETURN s, labels(s), s.cod_id, e, edgeType(e), c, labels(c), c.formula LIMIT 60`, @@ -464,6 +650,9 @@ RETURN s, labels(s), s.cod_id, e, edgeType(e), c, labels(c), c.formula LIMIT 60` }, { id: 'fragments', + shelf: 'solidform', + blurb: + 'Perceived functional groups are nodes, which is what makes synthon and coformer questions expressible at all.', prompt: 'Which functional groups dominate the corpus?', cypher: `MATCH (c:Component)-[e:HAS_FRAGMENT]->(f:Fragment) RETURN c, labels(c), c.formula, e, edgeType(e), f, labels(f), f.fragment_type LIMIT 400`, @@ -510,6 +699,9 @@ RETURN c, labels(c), c.formula, e, edgeType(e), f, labels(f), f.fragment_type LI }, { id: 'elements', + shelf: 'packing', + blurb: + 'Composition as traversal: element membership is an edge, so "organics containing both S and Br" needs no LIKE over a formula string.', prompt: 'What elements do these structures contain?', cypher: `MATCH (s:Structure)-[e:CONTAINS_ELEMENT]->(el:Element) RETURN s, labels(s), s.cod_id, e, edgeType(e), el, labels(el), el.symbol LIMIT 420`, @@ -535,6 +727,9 @@ RETURN s, labels(s), s.cod_id, e, edgeType(e), el, labels(el), el.symbol LIMIT 4 }, { id: 'solvates', + shelf: 'solidform', + blurb: + 'A hydrate is the same principal component with one extra edge -- so "every solvate of this compound" is one hop.', prompt: 'Show the solvates and hydrates', cypher: `MATCH (s:Structure)-[e:CONTAINS_COMPONENT]->(c:Component) WHERE c.is_solvent = true @@ -572,6 +767,9 @@ RETURN s, labels(s), s.cod_id, e, edgeType(e), c, labels(c), c.formula LIMIT 320 }, { id: 'acid_synthon', + shelf: 'solidform', + blurb: + 'The acid/carboxylate split distinguishes a co-crystal from a salt, which is a regulatory distinction rather than a cosmetic one.', prompt: 'Find the carboxylic acid molecules (synthon candidates)', cypher: `MATCH (c:Component)-[e:HAS_FRAGMENT]->(f:Fragment) WHERE f.fragment_type = 'carboxylic_acid' OR f.fragment_type = 'carboxylate' OR f.fragment_type = 'pyridine_nitrogen' @@ -603,13 +801,16 @@ RETURN c, labels(c), c.formula, e, edgeType(e), f, labels(f), f.fragment_type LI return [ `${acid.toLocaleString()} components carry a neutral carboxylic acid, ${carboxylate.toLocaleString()} carry the deprotonated carboxylate, and ${pyridine.toLocaleString()} carry a pyridine-type nitrogen acceptor. The canvas shows three clusters, one per fragment.`, `The acid/carboxylate split is doing real work. It distinguishes a co-crystal from a salt, which is a regulatory distinction in pharmaceutical solid-form work, not a cosmetic one -- and it falls out of the graph because proton position was perceived per component rather than assumed.`, - `From here the synthon question is a pattern match rather than a geometry job: take two acid-bearing components and ask for a reciprocal pair of hydrogen bonds between them, which is the R2,2(8) acid dimer. There are ${hbonds.toLocaleString()} hydrogen bonds with a located hydrogen available to satisfy it.`, - `Worth noting where the dialect bites: TuringDB rejects a cycle written as a cycle, so that reciprocal pair has to be expressed as two independent edge patterns joined on properties. The result is identical; the phrasing is not what a Neo4j author would write.`, + `From here the synthon question is a pattern match rather than a geometry job. There are ${hbonds.toLocaleString()} hydrogen bonds with a located hydrogen to work with; the "Find the carboxylic acid dimer motif" investigation takes it the rest of the way.`, + `One correction worth carrying, because it is the kind of mistake this representation invites. The obvious way to look for an R2,2(8) dimer is to ask for a RECIPROCAL PAIR of hydrogen bonds between two molecules. In this graph that pair does not exist: the stored network is a quotient graph, so the second bond of a centrosymmetric dimer is the symmetry image of the first and both collapse to ONE edge. Asking for a pair therefore matches each edge against itself. What actually separates a dimer from a chain is whether the generating operation is its own inverse -- and that is a boolean on the edge, not a second pattern.`, ] }, }, { id: 'history', + shelf: 'provenance', + blurb: + 'One commit per publication year: a statistic published in 2010 is a checkout, not an archived dump restore.', prompt: 'How has the corpus grown, and can I reproduce an old statistic?', graph: LEDGER_GRAPH, cypher: `MATCH (a:Snapshot)-[e:NEXT]->(b:Snapshot) @@ -617,21 +818,21 @@ RETURN a, labels(a), a.tag, e, edgeType(e), b, labels(b), b.tag`, answer: async () => { const L = LEDGER_GRAPH const V = VERSIONED_GRAPH - const [nSnap, firstYear, lastYear, totalNow] = await Promise.all([ + // This dialect implements only count() and avg() -- min, max and sum are + // all PLAN_ERROR -- so the ends of the range are read off the ordered + // chain rather than aggregated. + const [nSnap, totalNow] = await Promise.all([ scalar(L, 'MATCH (s:Snapshot) RETURN count(s)'), - scalar(L, 'MATCH (s:Snapshot) RETURN min(s.year)').catch(() => 0), - scalar(L, 'MATCH (s:Snapshot) RETURN max(s.year)').catch(() => 0), scalar(V, 'MATCH (s:Structure) RETURN count(s)'), ]) - // min/max do not exist in this dialect, so read the ends off the chain const ends = await runCypher( L, 'MATCH (s:Snapshot) RETURN s.year, s.tag, s.structures_total ORDER BY s.year' ) const years = (ends[0] ?? []).map((v) => Number(v)) const totals = (ends[2] ?? []).map((v) => Number(v)) - const y0 = years[0] ?? firstYear - const y1 = years[years.length - 1] ?? lastYear + const y0 = years[0] ?? 0 + const y1 = years[years.length - 1] ?? 0 const peakIdx = totals.reduce( (best, _v, i) => i > 0 && totals[i] - totals[i - 1] > totals[best] - totals[best - 1] ? i : best, @@ -647,6 +848,9 @@ RETURN a, labels(a), a.tag, e, edgeType(e), b, labels(b), b.tag`, }, { id: 'audit', + shelf: 'provenance', + blurb: + 'If a determination is later corrected, the commit history names exactly which published analyses included it.', prompt: 'If a structure is later corrected, which analyses are suspect?', graph: LEDGER_GRAPH, cypher: `MATCH (a:Snapshot)-[e:NEXT]->(b:Snapshot) @@ -679,10 +883,27 @@ RETURN a, labels(a), a.tag, e, edgeType(e), b, labels(b), b.tag`, }, { id: 'spacegroups', + shelf: 'packing', + blurb: + 'Symmetry as a node you can traverse to, not a text column -- which is why P2(1)/c and P2(1)/n stay distinguishable.', prompt: 'How are space groups distributed in this corpus?', cypher: `MATCH (s:Structure)-[e:IN_SPACE_GROUP]->(g:SpaceGroup) RETURN s, labels(s), s.cod_id, e, edgeType(e), g, labels(g), g.hm_symbol LIMIT 120`, timeline: { + answerAt: async (commit, tag) => { + const [structures, groups, p21c, p1bar, centro] = await Promise.all([ + scalarAt(VERSIONED_GRAPH, 'MATCH (s:Structure) RETURN count(s)', commit), + scalarAt(VERSIONED_GRAPH, 'MATCH (n:SpaceGroup) RETURN count(n)', commit), + scalarAt(VERSIONED_GRAPH, "MATCH (s:Structure)-[:IN_SPACE_GROUP]->(x:SpaceGroup) WHERE x.hm_symbol = 'P 1 21/c 1' RETURN count(s)", commit), + scalarAt(VERSIONED_GRAPH, "MATCH (s:Structure)-[:IN_SPACE_GROUP]->(x:SpaceGroup) WHERE x.hm_symbol = 'P -1' RETURN count(s)", commit), + scalarAt(VERSIONED_GRAPH, "MATCH (s:Structure)-[:IN_SPACE_GROUP]->(x:SpaceGroup) WHERE x.hm_symbol = 'C 1 2/c 1' RETURN count(s)", commit), + ]) + return [ + `At ${tag} the corpus held ${structures.toLocaleString()} structures over ${groups.toLocaleString()} Hermann-Mauguin settings.`, + `P 1 21/c 1 accounted for ${p21c.toLocaleString()} of them (${pct(p21c, structures)}%), P -1 for ${p1bar.toLocaleString()} (${pct(p1bar, structures)}%) and C 1 2/c 1 for ${centro.toLocaleString()} (${pct(centro, structures)}%).`, + `This is not a filtered view of today's corpus. It is the graph as it stood at that commit -- the structures published later do not exist in it. Drag along the strip and watch P2(1)/c's share climb: the skew that every crystallographer knows got MORE pronounced over the period, which is a claim you can only make if you can still run the query against the old corpus.`, + ] + }, metrics: async (commit) => { const [structures, groups, p21c, p1bar] = await Promise.all([ scalarAt(VERSIONED_GRAPH, 'MATCH (s:Structure) RETURN count(s)', commit), @@ -717,14 +938,231 @@ RETURN s, labels(s), s.cod_id, e, edgeType(e), g, labels(g), g.hm_symbol LIMIT 1 .filter((r) => r.n > 0) const nGroups = await scalar(g, 'MATCH (n:SpaceGroup) RETURN count(n)') const top = ranked.slice(0, 4).map((r) => `${r.hm} (${r.n.toLocaleString()}, ${pct(r.n, total)}%)`).join(', ') - const centro = ranked.filter((r) => ['P 1 21/c 1', 'P -1', 'C 1 2/c 1', 'P 1 21/n 1', 'P b c a'].includes(r.hm)) - .reduce((s2, r) => s2 + r.n, 0) + // Asked of the graph via the is_centrosymmetric flag, rather than summed + // from a hand-picked list of groups: the long tail of centrosymmetric + // groups is large enough that any short list understates it badly. + const centro = await scalar( + g, + 'MATCH (s:Structure)-[:IN_SPACE_GROUP]->(x:SpaceGroup) ' + + 'WHERE x.is_centrosymmetric = true RETURN count(s)' + ) return [ `${total.toLocaleString()} structures are distributed over ${nGroups} distinct space groups, and the distribution is extremely skewed.`, `The leaders are ${top}.`, - `Those centrosymmetric groups alone account for ${pct(centro, total)}% of the corpus. That is the familiar result: molecules pack most efficiently using inversion centres and glide planes, so a handful of groups dominates the whole of small-molecule crystallography.`, - `Note that P2(1)/c and P2(1)/n appear separately here. They are the same group (No. 14) in different cell settings, which is exactly the kind of distinction that gets lost when symmetry is a text column rather than a node you can traverse to.`, - `That distinction is also why two counts in this app differ, and it is worth being precise about. The panel on the left reports 159 space groups because that graph keys a SpaceGroup node by group NUMBER. The commit toggle reports more because the chronologically-committed graph keys by Hermann-Mauguin SYMBOL, so each setting is its own node. Same corpus, two modelling choices, and the node carries both properties so either question is answerable.`, + `Centrosymmetric groups account for ${pct(centro, total)}% of the corpus -- that figure is read off an is_centrosymmetric flag on the SpaceGroup node, not summed from the few groups on screen. It is the familiar result: molecules pack most efficiently using inversion centres and glide planes, so a handful of groups dominates the whole of small-molecule crystallography.`, + `One modelling point worth being exact about, because it changes what the numbers mean. This graph keys a SpaceGroup node by group NUMBER, so the ${(ranked[0]?.n ?? 0).toLocaleString()} structures shown under P2(1)/c are all of space group No. 14 in every setting -- P2(1)/n and P2(1)/a included. The chronologically-committed graph behind the commit toggle keys by Hermann-Mauguin SYMBOL instead, so there the settings are separate nodes and the count is higher. Same corpus, two deliberate modelling choices, and the node carries both properties so either question stays answerable.`, + ] + }, + }, + // ---- Packing & motifs ------------------------------------------------- + { + id: 'dimensionality', + shelf: 'packing', + blurb: + 'Whether a packing is a dimer, a chain, a sheet or a framework is the rank of its cycle lattice -- a traversal, with no relational form.', + prompt: 'Is this packing a dimer, a chain, a sheet or a framework?', + cypher: `MATCH (s:Structure)-[e:IN_SPACE_GROUP]->(g:SpaceGroup) +WHERE s.net_dim > 0 +RETURN s, labels(s), s.cod_id, e, edgeType(e), g, labels(g), g.hm_symbol LIMIT 140`, + answer: async (graph) => { + const g = graph + const [d0, d1, d2, d3, w0, w1, w2, w3, scored] = await Promise.all([ + scalar(g, 'MATCH (s:Structure) WHERE s.net_dim = 0 RETURN count(s)'), + scalar(g, 'MATCH (s:Structure) WHERE s.net_dim = 1 RETURN count(s)'), + scalar(g, 'MATCH (s:Structure) WHERE s.net_dim = 2 RETURN count(s)'), + scalar(g, 'MATCH (s:Structure) WHERE s.net_dim = 3 RETURN count(s)'), + scalar(g, 'MATCH (s:Structure) WHERE s.net_dim_weak = 0 RETURN count(s)'), + scalar(g, 'MATCH (s:Structure) WHERE s.net_dim_weak = 1 RETURN count(s)'), + scalar(g, 'MATCH (s:Structure) WHERE s.net_dim_weak = 2 RETURN count(s)'), + scalar(g, 'MATCH (s:Structure) WHERE s.net_dim_weak = 3 RETURN count(s)'), + scalar(g, 'MATCH (s:Structure) WHERE s.has_net_dim = true RETURN count(s)'), + ]) + const tot = d0 + d1 + d2 + d3 + const wtot = w0 + w1 + w2 + w3 + return [ + `Dimensionality is the rank of the lattice-translation subgroup generated by cycles in the hydrogen-bond net. Rank 0 is a finite motif -- a dimer or an isolated cluster. Rank 1 is a chain, rank 2 a sheet, rank 3 a framework.`, + `Over ${scored.toLocaleString()} structures carrying at least one located-hydrogen bond: ${d0.toLocaleString()} are 0D (${pct(d0, tot)}%), ${d1.toLocaleString()} 1D (${pct(d1, tot)}%), ${d2.toLocaleString()} 2D (${pct(d2, tot)}%) and ${d3.toLocaleString()} 3D (${pct(d3, tot)}%).`, + `Admitting weak C-H...A contacts as well moves it to ${pct(w0, wtot)}% 0D, ${pct(w1, wtot)}% 1D, ${pct(w2, wtot)}% 2D and ${pct(w3, wtot)}% 3D -- which SATURATES. Once every C-H donor counts, almost everything percolates and the statistic stops discriminating. That is the honest reading: the strong-bond net is the informative one, and the weak net is a reminder that a dimensionality number means nothing without the criterion attached to it.`, + `This matters because dimensionality predicts behaviour a formulator cares about: chains tend to needle morphology and anisotropic mechanical response, sheets to plates that cleave and tablet well, frameworks to harder, less soluble, humidity-stable solids.`, + `The mechanism is the reason this belongs in a graph, and it is also where the first version of this panel was wrong. The stored network is a QUOTIENT graph -- nodes are asymmetric-unit sites, each edge carries the symmetry operation to its neighbour -- so a chain running through the crystal is a short cycle, not a long path, and counting reachable nodes by depth reports everything as isolated. But the hydrogen bond alone is not the net: a donor and an acceptor are joined in the crystal THROUGH the molecules they belong to. Building the quotient graph from contacts only makes it bipartite with no path longer than one hop, so it has no cycles and everything comes back 0D. The covalent bonds have to be in the graph too, contracting each molecule to a point. A ring inside a finite molecule closes with the zero lattice vector, so it contributes nothing to the rank.`, + ] + }, + }, + { + id: 'acid_dimer', + shelf: 'packing', + blurb: + 'The R2,2(8) carboxylic-acid dimer, separated from the C(4) catemer by the ORDER of the symmetry operation that generates it.', + prompt: 'Find the carboxylic acid dimer motif, R2,2(8)', + cypher: `MATCH (f1:Fragment)<-[:IN_FRAGMENT]-(a1:Atom) + -[h:CONTACT]->(a2:Atom)-[:IN_FRAGMENT]->(f2:Fragment) +WHERE f1.fragment_type = 'carboxylic_acid' AND f2.fragment_type = 'carboxylic_acid' + AND h.kind = 'hbond' AND h.h_inferred = false + AND h.is_involution = true + AND a1.element = 'O' AND a2.element = 'O' + AND a1.cod_id = a2.cod_id +RETURN a1, labels(a1), a1.label, h, edgeType(h), a2, labels(a2), a2.label LIMIT 300`, + answer: async (graph, sub) => { + const g = graph + const [dimers, catemers, acids, acidStructures, meanLen, meanAng] = await Promise.all([ + scalar(g, `${DIMER_COUNT}true AND a1.element = 'O' AND a2.element = 'O' AND a1.cod_id = a2.cod_id RETURN count(h)`), + scalar(g, `${DIMER_COUNT}false AND a1.element = 'O' AND a2.element = 'O' AND a1.cod_id = a2.cod_id RETURN count(h)`), + scalar(g, "MATCH (c:Component)-[:HAS_FRAGMENT]->(f:Fragment) WHERE f.fragment_type = 'carboxylic_acid' RETURN count(c)"), + // Distinct molecular IDENTITIES and the structures they appear in are + // very different numbers here, and saying only the first invites the + // wrong reading. Components are deduplicated corpus-wide, and the + // common coformers repeat heavily: ~900 distinct acids across ~10,000 + // structures. That ratio IS the cross-structure identity argument. + scalar(g, "MATCH (f:Fragment)<-[:HAS_FRAGMENT]-(c:Component)<-[:CONTAINS_COMPONENT]-(s:Structure) WHERE f.fragment_type = 'carboxylic_acid' RETURN count(s)"), + scalar(g, `${DIMER_COUNT}true AND a1.element = 'O' AND a2.element = 'O' AND a1.cod_id = a2.cod_id RETURN avg(h.length)`), + scalar(g, `${DIMER_COUNT}true AND a1.element = 'O' AND a2.element = 'O' AND a1.cod_id = a2.cod_id RETURN avg(h.angle)`), + ]) + const tot = dimers + catemers + return [ + `${dimers.toLocaleString()} acid-to-acid hydrogen bonds are generated by an operation of order two, which is the R2,2(8) dimer. Mean H...O ${f2(meanLen)} A at ${f1(meanAng)} degrees -- textbook dimer geometry.`, + `Those come from ${acids.toLocaleString()} DISTINCT acid-bearing molecules appearing across ${acidStructures.toLocaleString()} structures, and every hit is an oxygen that is itself part of a carboxyl group -- the IN_FRAGMENT edge records which atoms are in the group, not merely which molecules contain one. The gap between those two numbers is the point: a Component node is deduplicated corpus-wide by identity, so the handful of coformers that crystal engineers actually reach for -- benzoic, maleic, succinic, fumaric -- are each one node carrying edges to every structure they appear in. That is what makes "where else has this molecule been co-crystallised" a traversal instead of a string join.`, + `The discriminator is the one piece of crystallography this whole demo rests on. A dimer is generated by an INVOLUTION -- an inversion centre, a mirror or a two-fold -- an operation that is its own inverse and therefore relates exactly two molecules. A catemer is generated by a 2(1) screw or a glide, which has infinite order and builds a chain instead. Both have a two-fold rotation part, so nothing short of composing the operation with itself tells them apart: ${catemers.toLocaleString()} acid-to-acid bonds here are catemer-generating and are excluded.`, + `That test is precomputed on every contact edge at ingest, so separating a ring motif from a chain motif costs one boolean in a WHERE clause rather than a symmetry analysis per candidate.`, + `Worth being precise about the representation, because it is the part people get wrong. The stored graph is a QUOTIENT graph: the second hydrogen bond of a centrosymmetric dimer is the symmetry image of the first, so the dimer is ONE edge here, not two. The R2,2(8) ring is a cycle in the unfolded net, closed by the operation the edge carries. The canvas shows ${sub.nodes.length} oxygen sites over ${sub.edges.length} such edges.`, + ] + }, + }, + { + id: 'weak_strong', + shelf: 'packing', + blurb: + 'Three contact populations kept deliberately separate: strong, heavy-atom-inferred, and weak C-H...A. Mixing them silently is how packing statistics go wrong.', + prompt: 'How do weak C-H...O contacts change the picture?', + cypher: `MATCH (a:Atom)-[e:CONTACT]->(b:Atom) +WHERE a.cod_id = ${PACKING_SEED_COD} +RETURN a, labels(a), a.label, e, edgeType(e), b, labels(b), b.label`, + answer: async (graph) => { + const g = graph + const [strong, inferred, weak, halogen, structures, mwLen, mwAng] = + await Promise.all([ + scalar(g, "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' AND r.h_inferred = false RETURN count(r)"), + scalar(g, "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond' AND r.h_inferred = true RETURN count(r)"), + scalar(g, "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond_weak' RETURN count(r)"), + scalar(g, "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'halogen' RETURN count(r)"), + scalar(g, 'MATCH (s:Structure) RETURN count(s)'), + scalar(g, "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond_weak' RETURN avg(r.length)"), + scalar(g, "MATCH ()-[r:CONTACT]->() WHERE r.kind = 'hbond_weak' RETURN avg(r.angle)"), + ]) + const total = strong + inferred + weak + halogen + return [ + `The contact layer holds ${total.toLocaleString()} edges over ${structures.toLocaleString()} structures, and they are three different claims, stored as three different populations.`, + `${strong.toLocaleString()} are strong hydrogen bonds with a located hydrogen. ${inferred.toLocaleString()} are flagged h_inferred -- no hydrogen was refined, so the criterion fell back to a heavy-atom donor-acceptor separation. ${weak.toLocaleString()} are weak C-H...A contacts, mean H...A ${f2(mwLen)} A at ${f1(mwAng)} degrees. ${halogen.toLocaleString()} are halogen bonds.`, + `The weak population is the one that changes conclusions. With only strong donors most packings do not percolate at all, which is why the dimensionality investigation reports so much 0D. C-H...O is a real, structure-directing interaction -- Taylor and Kennard established that from the CSD itself in 1982 -- and leaving it out makes a corpus look far more like isolated dimers than it is.`, + `The reason they are separate edge kinds rather than one blended "contact" is that a query must be able to refuse the weaker evidence. Every strong-bond statistic in this demo filters kind = 'hbond', so admitting the weak class changed none of them. A schema that merged them would have silently moved every number on the panel.`, + ] + }, + }, + // ---- Solid form & design ---------------------------------------------- + { + id: 'synthon_competition', + shelf: 'solidform', + blurb: + 'When a molecule offers two competing hydrogen-bond partners, which motif actually forms? A frequency count over the contact graph, conditioned on perceived fragments.', + prompt: 'When acid and pyridine compete, which synthon wins?', + cypher: `MATCH (f1:Fragment)<-[e:HAS_FRAGMENT]-(c:Component)-[:HAS_FRAGMENT]->(f2:Fragment) +WHERE f1.fragment_type = 'carboxylic_acid' AND f2.fragment_type = 'pyridine_nitrogen' +RETURN c, labels(c), c.formula, e, edgeType(e), f1, labels(f1), f1.fragment_type LIMIT 200`, + answer: async (graph, sub) => { + const g = graph + const F = (t: string) => + `MATCH (c:Component)-[:HAS_FRAGMENT]->(f:Fragment) WHERE f.fragment_type = '${t}' RETURN count(c)` + // `homo` is split by is_involution deliberately. A centrosymmetric acid + // homodimer is TWO physical hydrogen bonds stored as ONE quotient edge + // (the second is the symmetry image of the first), while an acid...N + // heterosynthon is one bond and one edge. Counting edges therefore + // understates the homosynthon by a factor of two, and on this corpus + // that is not academic: raw edge counts give 47.9% homo / 52.1% hetero, + // and correcting the double-counting REVERSES the winner. Rather than + // silently doubling, both figures are reported and the reason is stated. + const [both, acid, pyr, homoInv, homoPlain, hetero] = await Promise.all([ + scalar(g, "MATCH (f1:Fragment)<-[:HAS_FRAGMENT]-(c:Component)-[:HAS_FRAGMENT]->(f2:Fragment) WHERE f1.fragment_type = 'carboxylic_acid' AND f2.fragment_type = 'pyridine_nitrogen' RETURN count(c)"), + scalar(g, F('carboxylic_acid')), + scalar(g, F('pyridine_nitrogen')), + scalar(g, "MATCH (f1:Fragment)<-[:IN_FRAGMENT]-(a1:Atom)-[h:CONTACT]->(a2:Atom)-[:IN_FRAGMENT]->(f2:Fragment) WHERE f1.fragment_type = 'carboxylic_acid' AND f2.fragment_type = 'carboxylic_acid' AND h.kind = 'hbond' AND h.h_inferred = false AND h.is_involution = true AND a1.element = 'O' AND a2.element = 'O' RETURN count(h)"), + scalar(g, "MATCH (f1:Fragment)<-[:IN_FRAGMENT]-(a1:Atom)-[h:CONTACT]->(a2:Atom)-[:IN_FRAGMENT]->(f2:Fragment) WHERE f1.fragment_type = 'carboxylic_acid' AND f2.fragment_type = 'carboxylic_acid' AND h.kind = 'hbond' AND h.h_inferred = false AND h.is_involution = false AND a1.element = 'O' AND a2.element = 'O' RETURN count(h)"), + scalar(g, "MATCH (f1:Fragment)<-[:IN_FRAGMENT]-(a1:Atom)-[h:CONTACT]->(a2:Atom)-[:IN_FRAGMENT]->(f2:Fragment) WHERE f1.fragment_type = 'carboxylic_acid' AND f2.fragment_type = 'pyridine_nitrogen' AND h.kind = 'hbond' AND h.h_inferred = false AND a1.element = 'O' AND a2.element = 'N' RETURN count(h)"), + ]) + const homoEdges = homoInv + homoPlain + // an involution-generated edge stands for two physical hydrogen bonds + const homoBonds = homoInv * 2 + homoPlain + const totEdges = homoEdges + hetero + const totBonds = homoBonds + hetero + return [ + `${acid.toLocaleString()} components carry a carboxylic acid and ${pyr.toLocaleString()} a pyridine-type nitrogen. ${both.toLocaleString()} carry BOTH, and those are the interesting ones: the molecule offers the crystal a choice between bonding acid-to-acid and acid-to-pyridine.`, + `Counting stored edges, with inferred hydrogens excluded: ${homoEdges.toLocaleString()} acid-to-acid (${pct(homoEdges, totEdges)}%) against ${hetero.toLocaleString()} acid-to-pyridine (${pct(hetero, totEdges)}%). Counting physical HYDROGEN BONDS, it is ${homoBonds.toLocaleString()} (${pct(homoBonds, totBonds)}%) against ${hetero.toLocaleString()} (${pct(hetero, totBonds)}%).`, + `The two readings differ because ${homoInv.toLocaleString()} of the acid-to-acid edges are generated by an involution, and each of those stands for TWO physical hydrogen bonds: the second is the symmetry image of the first, so the quotient graph stores it once. The heterosynthon has no such pairing. Counting edges therefore understates the homodimer, which is the kind of bias a stored-contact representation introduces and that has to be corrected out loud rather than quietly. Here it widens the margin rather than changing the answer -- the homodimer wins on either reading.`, + `One scoping caveat worth saying aloud: only ${both.toLocaleString()} molecules in this corpus carry BOTH groups, so these are corpus-wide propensities, not a head-to-head contest within molecules that had a real choice.`, + `This is the question co-crystal design turns on. The acid/pyridine heterosynthon is the workhorse of pharmaceutical co-crystal screening precisely because it tends to beat the acid homodimer, and a corpus-wide frequency is the evidence a formulator would want before committing screening time.`, + `The canvas shows ${sub.nodes.length} nodes. Note what had to exist for this to be one query: perceived fragments as NODES, an atom-to-molecule edge so a contact between two atoms resolves to a contact between two molecules, and the hydrogen-bond population stored rather than recomputed. Take away any one and this becomes a scripted job over a structure archive.`, + ] + }, + }, + { + id: 'polymorph', + shelf: 'solidform', + blurb: + 'The same molecule packing two different ways. Only expressible because molecular identity is one node shared across structures.', + prompt: 'Which molecules crystallise in more than one packing?', + cypher: `MATCH (s:Structure)-[e:CONTAINS_COMPONENT]->(c:Component) +WHERE e.role = 'principal' AND c.is_solvent = false AND c.has_inchikey = true +RETURN c, labels(c), c.formula, e, edgeType(e), s, labels(s), s.hm_symbol LIMIT 260`, + answer: async (graph, sub) => { + const g = graph + // Keyed on space-group NUMBER, not Hermann-Mauguin symbol. An earlier + // version compared hm_symbol, which counted P2(1)/c against P2(1)/n -- + // the same group in a different setting -- as a different packing, and + // so contradicted the space-group panel two clicks away. It also counted + // components whose identity is a connectivity digest rather than an + // InChIKey; both are now excluded. + // Written as one connected chain through the SpaceGroup nodes, which is + // both how the group NUMBER is reachable and, in this engine, ~20x + // faster than joining the two structures as separate comma patterns. + const base = + "MATCH (g1:SpaceGroup)<-[:IN_SPACE_GROUP]-(s1:Structure)" + + "-[r1:CONTAINS_COMPONENT]->(c:Component)" + + "<-[r2:CONTAINS_COMPONENT]-(s2:Structure)-[:IN_SPACE_GROUP]->(g2:SpaceGroup) " + + "WHERE r1.role = 'principal' AND r2.role = 'principal' " + + "AND c.is_solvent = false AND c.has_inchikey = true " + const [pairs, sameGroup, comps, withKey] = await Promise.all([ + scalar(g, `${base}AND g1.number <> g2.number RETURN count(c)`).catch(() => 0), + scalar(g, `${base}AND g1.number = g2.number AND s1.cod_id <> s2.cod_id RETURN count(c)`).catch(() => 0), + scalar(g, 'MATCH (c:Component) RETURN count(c)'), + scalar(g, 'MATCH (c:Component) WHERE c.has_inchikey = true RETURN count(c)'), + ]) + return [ + `${pairs.toLocaleString()} ordered pairs of structures share a principal molecule of identical InChIKey but crystallise in a different space GROUP -- not merely a different setting of the same group, which is the distinction the space-group panel draws.`, + `A further ${sameGroup.toLocaleString()} pairs share both the molecule and the space group. Those are where the interesting cases hide: two forms in the same group with different cells are still polymorphs, and no space-group comparison will ever find them. Separating the two counts is the honest way to present it.`, + `The reason either is a one-hop question is the modelling choice underneath. A Component node is deduplicated corpus-wide by identity, so the same molecule appearing in forty structures is ONE node with forty CONTAINS_COMPONENT edges. In a relational schema this is a self-join over a structure table on a chemical identity string -- which only works if that string was computed consistently in the first place, and still says nothing about the packing without a second pass.`, + `Two limits stated rather than buried. Identity is an InChIKey where RDKit perception succeeds, ${pct(withKey, comps)}% of ${comps.toLocaleString()} components; metal complexes fail by construction, since InChI has no well-defined representation for them, and those carry a connectivity digest in a SEPARATE property so it is never mistaken for a key. And nothing here requires the rest of the composition to match, so a hydrate against an anhydrate is included -- strictly those are solvates, and the distinction is regulatory. The canvas shows ${sub.nodes.length} nodes.`, + ] + }, + }, + { + id: 'coformer', + shelf: 'solidform', + blurb: + 'Molecules that co-crystallise with the partners your target already has, but have never been tried with the target itself -- a five-hop join plus an anti-join.', + prompt: 'Suggest coformers that have never been tried with this compound', + cypher: `MATCH (f:Fragment)<-[:HAS_FRAGMENT]-(c1:Component)<-[r1:CONTAINS_COMPONENT]-(s:Structure)-[e:CONTAINS_COMPONENT]->(c2:Component) +WHERE f.fragment_type = 'carboxylic_acid' AND c2.is_solvent = false +RETURN c2, labels(c2), c2.formula, e, edgeType(e), s, labels(s), s.cod_id LIMIT 240`, + answer: async (graph, sub) => { + const g = graph + const [partners, multi, solvates] = await Promise.all([ + scalar(g, "MATCH (f:Fragment)<-[:HAS_FRAGMENT]-(c1:Component)<-[:CONTAINS_COMPONENT]-(s:Structure)-[:CONTAINS_COMPONENT]->(c2:Component) WHERE f.fragment_type = 'carboxylic_acid' AND c2.is_solvent = false RETURN count(c2)"), + scalar(g, 'MATCH (s:Structure) WHERE s.n_components > 1 RETURN count(s)'), + scalar(g, "MATCH (s:Structure)-[r:CONTAINS_COMPONENT]->(c:Component) WHERE c.is_solvent = true RETURN count(r)"), + ]) + return [ + `${partners.toLocaleString()} non-solvent component pairings sit in structures where one partner carries a carboxylic acid. ${multi.toLocaleString()} structures are multi-component, and ${solvates.toLocaleString()} component slots are occupied by a classified solvent, which is why the solvent flag has to be in the query rather than applied afterwards -- otherwise every suggestion comes back as water.`, + `The commercial shape of this question is: my compound has these hydrogen-bond donors and acceptors. Find molecules that are known to co-crystallise with compounds carrying the same complement, but that have never been tried with mine. That is five hops out -- target to fragment, fragment to other molecules carrying it, those to their structures, structures to their partners -- followed by an ANTI-join to remove the pairs already reported.`, + `The traversal is the part a graph does well and a relational schema does not; the anti-join and the ranking are done client-side here because this dialect has no collect() or WITH, which is a real limitation worth naming rather than hiding.`, + `The canvas shows ${sub.nodes.length} nodes over ${sub.edges.length} co-occurrence edges. What is on screen is the evidence layer -- who has actually been crystallised with whom -- which is the part nobody currently stores as a graph.`, ] }, },