diff --git a/utilities/diagram_gen/README.md b/utilities/diagram_gen/README.md new file mode 100644 index 0000000000..962a466338 --- /dev/null +++ b/utilities/diagram_gen/README.md @@ -0,0 +1,282 @@ +# diagram_gen — Goldstone diagrams from SeQuant term strings + +Turn a SeQuant coupled-cluster term into a publication-style Goldstone diagram (PNG). + +``` +python3 csv_diagram.py '' out.png +``` + +Everything is driven by the **term string** that SeQuant prints — you paste it in, you +get a diagram out. The renderer understands the CSV/PNO–DF flavour of those strings +(PAO→PNO transforms `C`, PAO overlaps `s`, density-fitting auxiliary index `Κ`) and can +either draw them as annotations or drop them to leave a plain textbook CC diagram. + +--- + +## 1. Requirements + +Python 3 with `matplotlib` and `numpy`. Nothing else; no SeQuant build, no LaTeX +installation (matplotlib's built-in mathtext renders the labels). + +```bash +python3 -c "import matplotlib, numpy; print('ok')" +``` + +--- + +## 2. Quick start + +```bash +cd utilities/diagram_gen + +# a) one term, given inline +python3 csv_diagram.py '2 (f{i_1;μ̃_1} * C{μ̃_1;a_1}) * t{a_1;i_1}' bubble.png + +# b) one term, from a file +python3 csv_diagram.py --file term.txt out.png + +# c) many terms: one file, one PNG per line -> term001.png, term002.png, ... +python3 csv_diagram.py --file sample_terms.txt outdir/ + +# d) the canonical CC skeleton instead of the CSV/DF overlay +python3 csv_diagram.py --file term.txt out.png --no-csv +``` + +`sample_terms.txt` in this folder holds four real terms to try; `examples/` shows what +they render to. + +| example | what it shows | +|---|---| +| `examples/ccsd_E3_fock_bubble.png` | one-body `—×` vertex, T1, particle+hole drawn as a lens | +| `examples/ccsd_R2t49_vvvv_ladder.png` | (vv\|vv) ladder + T2, four `C` boxes, `K` on the interaction | +| `examples/ccsd_R2t30_oooo_csc_bridge.png` | (oo\|oo) + two T1, `CsC` domain-changing bridges | +| `examples/ccsdt_R1t1_with_T3.png` | a T3 amplitude bar (three excitation apexes) | +| `examples/ccsd_R2t49_canonical.png` | same term as the ladder above, with `--no-csv` | + +--- + +## 3. The input: what a term string looks like + +A term is a product of tensors, written `label{bra;ket}` (or `label{bra;ket;aux}`), +joined by `*`. Parentheses are irrelevant — the renderer ignores the association. + +``` +1/2 Ŝ{i_1,i_2;a_1,a_2} * g{μ̃_1;μ̃_2;Κ_1} * C{a_1;μ̃_1} * … * t{a_3,a_4;i_1,i_2} + ↑ coefficient ↑ external indices (projector) ↑ interaction ↑ PAO→PNO transform ↑ amplitude +``` + +### Tensors understood + +| label | slots | meaning | drawn as | +|---|---|---|---| +| `g{p;q;Κ}` | bra, ket, auxiliary | one half of a two-body interaction (a DF three-centre factor) | a dot; the two `g`'s sharing a `Κ` are joined by the dashed interaction line | +| `f{p;q}` | bra, ket | one-body (Fock) operator | a dot joined by a dashed line to a `×` | +| `t{virt…;occ…}` | bra=virtuals, ket=occupieds | cluster amplitude T₁…T₄ | solid horizontal bar at the bottom, one V-shaped particle/hole pair per excitation | +| `C{p;q}` | bra, ket | PAO→PNO transform | *not* a vertex — it splices two indices into one fermion line and puts a `C` box on it | +| `s{μ̃;μ̃}` | bra, ket | PAO–PAO overlap metric | same, part of a `CsC` bridge | + +### Index naming (this is how line types are decided) + +| starts with | kind | drawn as | +|---|---|---| +| `i…` | occupied | **hole** line — blue, arrow pointing down | +| `a…` | PNO virtual | **particle** line — red, arrow pointing up | +| `μ̃…` | PAO virtual | internal to a `C`/`s` chain; never a line end in the picture | +| `Κ…` | DF auxiliary | the `K` label on the dashed interaction line | + +`a_3` means "virtual `a_3` restricted to pair domain (i₁i₂)"; the `<…>` part is +printed as the superscript `a₃^{i₁i₂}`. + +### What gets cleaned up automatically + +all of these are stripped for you: + +* `:N-C-S` (and similar) suffixes after a tensor +* markdown list numbering and backticks: ``12. `` `` +* surrounding whitespace, a leading `+` +* any parenthesisation + +So this works as a single argument, unchanged: + +``` +1. `-6 Ŝ{i_1,i_2,i_3;a_1,…}:N-C-S * g{i_4;μ̃_1;Κ_1}:N-C-S * …` +``` + +### Rules the term must satisfy + +1. **Every index appears exactly twice** — once in a bra slot, once in a ket slot. + A truncated paste is the usual cause of failure here. +2. **At most one interaction**: either two `g`'s sharing one `Κ`, or one `f`. + True of every CC residual term SeQuant emits. +3. Amplitudes may be T₁…T₄ (verified up to T₄; the bar just grows an apex per excitation). + +--- + +## 4. `csv_diagram.py` — render one term (or a small batch) + +``` +python3 csv_diagram.py [TERM] [OUT.png] [options] +python3 csv_diagram.py --file PATH [OUT.png|OUTDIR/] [options] +``` + +| option | effect | +|---|---| +| `TERM` | the term string; use `-` to read it from stdin | +| `OUT.png` | output path (default `diagram.png`). With `--file` and several terms, this is a **directory** and files are named `term001.png`, … | +| `--file PATH` | read term(s) from a file (`-` = stdin). Every line containing `{…}` is rendered; other lines (comments, blanks) are skipped | +| `--title TEXT` | figure title. The coefficient is appended automatically | +| `--no-csv` | **canonical CC skeleton**: drop the `C`/`CsC`/`s` boxes, the domain superscripts and the `K` / `g_L,g_R` tags | +| `--no-dom` | drop the domain superscripts only (`a₃^{i₁i₂}` → `a₃`) | +| `--no-aux` | drop the `K` label and the `g_L`/`g_R` tags only | +| `--no-legend` | omit the legend strip at the bottom | + +With no arguments at all it renders three built-in demo terms into the current directory. + +### The overlay knob + +`--no-csv` implements rule **B0** of `csv-cc-diagram-rules.md` (skeleton invariance): +deleting every domain superscript, `C`/`CsC` box and the `K` label must leave a textbook +Goldstone diagram. Compare `examples/ccsd_R2t49_vvvv_ladder.png` with +`examples/ccsd_R2t49_canonical.png` — same skeleton, no annotations. The three layers +(`--no-csv`, `--no-dom`, `--no-aux`) can be mixed freely. + +--- + +## 5. `render_all.py` — batch a whole equation file + +Renders every numbered term of an equation markdown file (sections `E`, `R1`, `R2`, +`R3`, `R4`) and writes an `INDEX.md` manifest next to the PNGs. + +```bash +# whole CCSD file, CSV overlay + canonical skeleton +python3 render_all.py --src pao-ccsd-df.md --variant both + +# just the triples section of the CCSDT file +python3 render_all.py --src pao-ccsdt-df.md --only R3 + +# three specific R4 terms into a scratch directory +python3 render_all.py --src pao-ccsdtq-df.md --only R4:1,2,3 --outdir /tmp/q +``` + +| option | effect | +|---|---| +| `--src FILE` | the equation markdown file. Defaults to `pao-ccsd-df.md` next to the script or one directory up | +| `--variant csv\|canonical\|both` | overlay, plain skeleton, or both (default `csv`) | +| `--only TAG[:N,N,…]` | restrict to a section (`R3`) or specific terms (`R2:16,49`) | +| `--outdir DIR` | output directory (single-variant runs). Default: `ccsd_all`, `ccsdt_all`, `ccsdtq_all`, … derived from the file name, with `_canonical` appended for the skeleton variant | +| `--kind TEXT` | title prefix (default derived from the file name, e.g. `CSV-CCSDT`) | + +It expects terms formatted as a markdown numbered list under a `## E`/`## R1`/… heading: + +```markdown +## R2 — doubles residual + +1. `` +2. `` +``` + +Failures never abort the run — the offending term is listed in `INDEX.md` with its error. + +--- + +## 6. `audit_diagrams.py` — check the output is clean + +Re-runs the layout for every term of a file and reports anything unreadable: overlapping +labels, labels sitting on a vertex or a bar, labels on a foreign line, and fermion lines +drawn on top of each other. + +```bash +python3 audit_diagrams.py --src pao-ccsdtq-df.md +python3 audit_diagrams.py --src pao-ccsd-df.md --no-csv +``` + +Over the CCSD/CCSDT/CCSDTQ equation files the renderer was developed against (660 terms, +both variants) the audit reports **no label/label, label/vertex or line/line overlaps**. +The only remaining reports are labels sitting on top of one *foreign* line — 1 of 86 +CCSD terms, 8 of 201 CCSDT, 32 of 373 CCSDTQ. Those labels keep an opaque background, +so they stay readable. + +--- + +## 7. Using it from Python + +```python +import csv_diagram as cd + +cd.draw(term_string, title="R2 term 49", out="out.png", + show_csv=True, # C / CsC / s boxes + show_dom=None, # domain superscripts (None = follow show_csv) + show_aux=None, # K and g_L/g_R tags (None = follow show_csv) + legend=True, dpi=170, verbose=False) # -> returns the output path +``` + +The parsing/graph layer is reusable on its own: + +```python +coeff, tensors = cd.parse_term(term_string) +nodes, lines, hvs, pair, ext = cd.build_graph(tensors) +for L in lines: + print("hole" if L.is_hole else "particle", + L.src, L.src_idx[0], "->", L.dst, L.dst_idx[0], L.annotation) +``` + +`build_graph` collapses every `C`/`s` chain, so `lines` is the physical fermion-line +graph: each `FermionLine` runs from the vertex where its index sits in a **bra** slot to +the one where it sits in a **ket** slot, and `.annotation` is `None`, `"C"` or `"CsC"` +depending on what the chain contained. + +--- + +## 8. How to read the diagram + +| element | meaning | +|---|---| +| red line, arrow up | particle line (virtual index) | +| blue line, arrow down | hole line (occupied index) | +| solid black bar | a cluster amplitude T_m; each excitation is one V of a particle + a hole leg meeting at an apex on the bar. The label under the bar is `t̂` with its occupied indices as subscripts and virtuals as superscripts | +| two grey dots joined by a dashed line | **one** two-body interaction vertex (its two half-vertices); `g_L`/`g_R` tag the left/right half | +| `K` on the dashed line | the density-fitting auxiliary index shared by the two `g` factors | +| dot + dashed line + `×` | the one-body (Fock) vertex | +| `C` box on a line | a PAO→PNO transform on that interaction leg | +| `CsC` box on a line | a domain-changing bridge C·s·C joining two PNOs of different domains | +| two lines bowed into a slim pointy ellipse | a particle and a hole line connecting the *same* two points — the canonical Goldstone bubble (they would otherwise be drawn on top of each other) | +| faint dotted line across the top | the projection boundary; lines crossing it are the external indices of `Ŝ` | + +Conventions and their textbook citations are in `csv-cc-diagram-rules.md`. Sign/weight +rules are **not** applied — the coefficient printed in the title is whatever the term +carried. + +--- + +## 9. When it fails + +Errors are raised as `ValueError` with a specific message: + +| message | cause | +|---|---| +| `no tensor found in input: …` | the string had no `{…}` — wrong argument, or the shell ate it (quote with single quotes) | +| `chain has N terminals (want 2)` | an index appears once, or three or more times — usually a truncated paste | +| `chain : both terminals are bra slots` | the same index used twice on the same side | +| `unknown index kind: 'x_1'` | an index name that starts with something other than `i`, `a`, `μ̃`, `Κ` | +| `>2 interaction half-vertices not supported` | more than one interaction in the term | + +Shell tip: always wrap the term in **single** quotes — it contains `*`, `{`, `}` and +backticks, all of which the shell would otherwise interpret. If quoting gets awkward, +put the term in a file and use `--file`. + +--- + +## 10. Files + +| file | role | +|---|---| +| `csv_diagram.py` | the renderer: parser, fermion-line graph, layout, drawing, CLI. Self-contained | +| `render_all.py` | batch driver over an equation markdown file; writes `INDEX.md` | +| `audit_diagrams.py` | machine check for overlapping labels/lines (imports the other two) | +| `csv-cc-diagram-rules.md` | the diagram conventions with citations (Shavitt & Bartlett, Crawford & Schaefer) | +| `sample_terms.txt` | four real terms to try | +| `examples/` | what those terms render to | + +`render_all.py` and `audit_diagrams.py` operate on equation markdown files +(`pao-ccsd-df.md`, `pao-ccsdt-df.md`, `pao-ccsdtq-df.md`) that are not part of this +repository; point them at your own with `--src`. diff --git a/utilities/diagram_gen/audit_diagrams.py b/utilities/diagram_gen/audit_diagrams.py new file mode 100644 index 0000000000..17224a8d10 --- /dev/null +++ b/utilities/diagram_gen/audit_diagrams.py @@ -0,0 +1,134 @@ +#!/usr/bin/env python3 +""" +audit_diagrams.py — machine check that every rendered CSV-CC diagram is clean. + +For each term of pao-ccsd-df.md it re-runs the layout and reports + 1. label/label overlaps (two text boxes sharing area) + 2. label/vertex-or-bar overlaps + 3. labels sitting on top of a fermion line that is not their own + 4. fermion lines drawn on top of each other (identical endpoints not bowed + apart, or two segments nearly collinear over a long stretch) + +Exit status 0 iff nothing is reported. Run for both variants: + python3 audit_diagrams.py # CSV overlay + python3 audit_diagrams.py --no-csv # canonical skeleton +""" + +import argparse +import os +import sys +import tempfile + +import numpy as np + +import csv_diagram as cd +from render_all import SRC, extract_terms, filter_sections, section_order + + +_TMP = os.path.join(tempfile.gettempdir(), "_csvdiag_audit.png") + + +def overlaps(a, b): + dx = min(a[2], b[2]) - max(a[0], b[0]) + dy = min(a[3], b[3]) - max(a[1], b[1]) + return dx * dy if (dx > 0 and dy > 0) else 0.0 + + +def audit_term(expr, show_csv): + cd.draw(expr, title=None, out=_TMP, verbose=False, + show_csv=show_csv, legend=False) + A = cd.LAST_AUDIT + boxes, tags, segs = A["boxes"], A["tags"], A["segs"] + issues = [] + + lab = [(b, t, ow) for b, (t, is_lab, ow) in zip(boxes, tags) if is_lab] + obs = [(b, t) for b, (t, is_lab, ow) in zip(boxes, tags) if not is_lab] + + for i in range(len(lab)): + for j in range(i + 1, len(lab)): + ov = overlaps(lab[i][0], lab[j][0]) + if ov > 1e-4: + issues.append(f"label/label {lab[i][1]} ~ {lab[j][1]} ({ov:.3f})") + for b, t, _ in lab: + for ob, ot in obs: + ov = overlaps(b, ob) + if ov > 1e-4: + issues.append(f"label/{ot} {t} ({ov:.3f})") + + # a label box sitting on a fermion line OTHER than the one it belongs to + for b, t, own in lab: + foreign = set() + for p, q, o in segs: + if o in ("dash", None) or id(o) == own: + continue + if cd._seg_hits_box(np.array(p), np.array(q), b): + foreign.add(id(o)) + if foreign: + issues.append(f"label sits on {len(foreign)} foreign line(s): {t}") + + # duplicated / overlapping fermion strokes: sample each line's polyline and + # measure how much of it runs within HAIR of another line + paths = {} + for p, q, o in segs: + if o in ("dash", None): + continue + paths.setdefault(id(o), []).append((np.array(p), np.array(q))) + sampled = {} + for k, ss in paths.items(): + pts = [ss[0][0]] + [b for _, b in ss] + pts = np.array(pts) + d = np.cumsum(np.r_[0, np.hypot(*np.diff(pts, axis=0).T)]) + if d[-1] < 1e-9: + continue + u = np.linspace(0, d[-1], 60) + sampled[k] = np.column_stack([np.interp(u, d, pts[:, 0]), + np.interp(u, d, pts[:, 1])]) + HAIR = 0.07 + keys = list(sampled) + for i in range(len(keys)): + for j in range(i + 1, len(keys)): + A, B = sampled[keys[i]], sampled[keys[j]] + dist = np.hypot(A[:, None, 0] - B[None, :, 0], A[:, None, 1] - B[None, :, 1]) + frac = float((dist.min(axis=1) < HAIR).mean()) + if frac > 0.35: + issues.append(f"two fermion lines coincide over {frac:.0%} of a line") + return issues + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--src", default=SRC, help="SeQuant equation markdown file") + ap.add_argument("--only", default=None, help="e.g. R3 or R4:1,2") + ap.add_argument("--no-csv", action="store_true") + args = ap.parse_args() + show_csv = not args.no_csv + + src = os.path.abspath(os.path.expanduser(args.src)) + if not os.path.exists(src): + ap.error(f"equation file not found: {src} (pass one with --src )") + sections = extract_terms(src) + if args.only: + sections = filter_sections(sections, args.only) + bad = n = 0 + for tag in section_order(sections): + for num, expr in sections.get(tag, []): + n += 1 + try: + issues = audit_term(expr, show_csv) + except Exception as exc: # noqa: BLE001 + print(f"{tag} term {num}: RENDER FAILED {type(exc).__name__}: {exc}") + bad += 1 + continue + if issues: + bad += 1 + print(f"{tag} term {num}:") + for it in issues: + print(f" {it}") + variant = "CSV overlay" if show_csv else "canonical" + print(f"\n[{variant}] {os.path.basename(args.src)}: " + f"{bad} of {n} term(s) with issues") + return 1 if bad else 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/utilities/diagram_gen/csv-cc-diagram-rules.md b/utilities/diagram_gen/csv-cc-diagram-rules.md new file mode 100644 index 0000000000..92decc4820 --- /dev/null +++ b/utilities/diagram_gen/csv-cc-diagram-rules.md @@ -0,0 +1,214 @@ +# CSV-CC Goldstone diagram rules + +Rule set for drawing and interpreting CSV-CCSD (cluster-specific-virtual, locally +density-fitted) Goldstone diagrams. Part A is the canonical textbook layer; Part B is +exactly what the CSV/DF overlay adds. Strip every Part-B annotation and a valid +CSV-CC diagram must reduce to the textbook diagram of Part A. + +Citations: +- **C&S** = Crawford & Schaefer 2000, *An Introduction to Coupled Cluster Theory*, + chapter "An Introduction to Coupled Cluster Diagrams", printed pp. 77–95 + (PDF page = printed − 32; formalism credited to Kucharski & Bartlett, p. 77). +- **S&B** = Shavitt & Bartlett 2009, *Many-Body Methods in Chemistry and Physics*, + Ch. 4 (printed page = PDF − 16), plus the CC rules of Ch. 9–10 (Fig. 10.1, p. 296). +- **KST** = K. S. Thanthiriwatte, *Diagrammatic techniques in Coupled-Cluster + Theory*, slides, Aug 2010 (18 slides). + +--- + +## Part A — canonical skeleton rules + +### A1. Line direction and labels +Upward-directed lines are **particle** lines (virtuals a, b, c, …); downward-directed +lines are **hole** lines (occupieds i, j, k, …). "Downward-directed lines represent +hole states (orbitals occupied in the reference) and upward directed lines represent +particle states" (C&S p. 78, Fig. 1; S&B p. 91: "pointing upward for particles and +downward for holes"; KST slide 4). The horizontal position of a line "has no +significance" (S&B p. 91). The reference |0⟩ is drawn as nothing (C&S Fig. 1c; +KST slide 4). + +### A2. Vertex inventory +- **One-body operator (Fock)**: a dot joined by a horizontal *dashed* line to a + marker "×" (C&S Fig. 2, p. 79; S&B p. 93 "the point of action … marked by the + interaction line —×"; KST slide 5). Four fragments f_ab, f_ij, f_ia, f_ai with + excitation levels 0, 0, −1, +1. +- **Two-body operator**: ONE vertex = **two half-vertices at the same level joined by + a horizontal dashed interaction line** ("The two half-vertices and the interaction + line constitute a single vertex", S&B p. 111, Eq. 4.28; C&S Fig. 3, p. 80: the nine + V_N fragments; KST slide 6). Electron assignment: "electron 1 ↔ left half-vertex; + electron 2 ↔ right half-vertex" (S&B p. 111, boxed). +- **Cluster operators T_m**: a **solid horizontal bar** with m hole–particle line + pairs rising from it; no lines below the bar ("the cluster operators contain only + q-creation strings … they contain no lines below the horizontal bar", C&S p. 80, + Fig. 4; S&B pp. 273–274, Eqs. 9.84–9.87 "solid horizontal lines to represent the + T̂_m operators"; KST slides 4, 8). + +### A3. Half-vertex degree rule +"**Each individual half-vertex will have one incoming and one outgoing line**, each of +which may be a particle line or a hole line" (S&B p. 111, verbatim). This holds for +every interaction dot (one-body or half of two-body). + +### A4. Bra = outgoing = creation; ket = incoming = annihilation +"incoming line ↔ annihilation operator ↔ ket state (IAK); outgoing line ↔ creation +operator ↔ bra state (OCB)" (S&B p. 98, mnemonics; also p. 95 rule 1 and the boxed +table on p. 111). Creation operators lie above the interaction line, annihilation +below (C&S p. 78; KST slide 4). Consequently: a particle line above a vertex with +arrow leaving it = created particle; a hole line below a vertex with arrow leaving +it = created hole; etc. + +### A5. Integral assignment +- One-body: **f = ⟨out | f | in⟩** — "(out|f|in), where out indicates the index of the + outgoing directed line and in the incoming" (C&S p. 84; KST slide 14 rule 2; + S&B p. 93: bra = leaving line, ket = entering line). +- Two-body: **⟨left-out right-out ‖ left-in right-in⟩**, with operator product + {(left-out)† (right-out)† (right-in)(left-in)} (S&B p. 112, boxed, verbatim; + C&S p. 85; KST slide 6, verbatim). Goldstone diagrams proper use the + non-antisymmetrized ⟨pq|rs⟩ (S&B p. 112); antisymmetrized-Goldstone (Brandow/ASG) + and Hugenholtz variants use ⟨pq‖rs⟩ (S&B pp. 118–123). *The CSV-CCSD spin-traced + equations are in the non-antisymmetrized (Goldstone-proper) regime.* + +### A6. Amplitude assignment +"With every T_m vertex associate an amplitude t^{ab…}_{ij…}" where the (hole, +particle) pairs are read off the bar **left to right**, each vertical pair supplying +one (subscript, superscript) column (KST slide 14 rule 4; C&S pp. 84–85 "hole and +particle indices in their left-to-right order in the diagram"; S&B p. 282 and +Fig. 10.1 rule 4). External pairing matches the bra determinant: "surviving hole and +particle lines at the top are paired vertically in |Φ…⟩ if they are part of the same +path" (S&B p. 104). + +### A7. Vertical placement (fictitious time) +Operators act bottom-to-top: rightmost operator (a T) at the bottom, Hamiltonian +fragment above it, resulting/bra state at the top (C&S p. 81, p. 83 "the rightmost +operator's interaction line must lie at the bottom"; S&B p. 91; KST slide 9). +T1 vertices commute, so the relative order/side-by-side placement of T bars is +immaterial (C&S p. 86; S&B p. 274, Eq. 9.88). External lines run to the top edge; +energy diagrams "can contain no external lines" (C&S p. 83); the R1/R2 equations have +exactly 2/4 external lines (C&S p. 88). + +### A8. Internal summation +"Summations are included over all 'internal' indices — lines that begin and end at +operator interaction lines and do not extend to infinity above or below" (C&S p. 84; +KST slide 14 rule 5: "Sum over all internal line labels, i.e. lines terminating +below the last Ĥ_N"; S&B p. 100, Eq. 4.19). + +### A9. Sign rule: (−1)^{h+l} with quasi-loops +Sign = (−1)^{h+l}, h = number of hole lines, l = number of loops. "A loop is a route +along a series of directed lines that either returns to its beginning or begins at +one external line and ends at another" (C&S p. 84). Open diagrams: "for open diagrams +fictitious external loop should be included" (KST slide 14 rule 6); S&B formalizes +this as **quasiloops** — "open lines … paired at the top according to the pairing +pattern in the corresponding Slater determinant" (S&B p. 135, Eq. 5.15; Fig. 10.1 +rule 8, p. 296). S&B writes the equivalent form (−1)^{h−l} (p. 109, Eq. 4.26). + +### A10. Equivalent lines → factor 1/2 +"A pair of lines is equivalent if they connect the same pair of vertices in the same +direction. Each pair … contributes a factor 1/2" (S&B p. 120; C&S p. 85: lines +"beginning at the same operator interaction line and ending at the same interaction +line"; KST slide 15 rule 7: weight (1/2)^m). Groups of n ≥ 3 equivalent lines +(MBPT/UCC only) give 1/n! (C&S p. 85, footnote q). + +### A11. Equivalent T vertices → factor 1/n! +"If there are n equivalent vertices in the diagram, they contribute a prefactor of +1/n!" (C&S p. 87). Two T vertices are equivalent when they have the same m and +"are connected to the interaction vertex with the same number … of particle lines +and the same number of hole lines" (S&B p. 287; KST slide 15 rule 9). + +### A12. Permutation operators for external lines +"Each pair of unique, external hole or particle lines introduces a permutation +function P(pq)", P(pq) f(p,q) = f(p,q) − f(q,p); pairs originating from the *same* +interaction line are not unique (already antisymmetric) (C&S p. 91, Eq. 154 p. 76; +KST slide 15 rule 8; S&B Fig. 10.1 rule 10). *In the spin-traced biorthogonal +CSV equations this role is played by the projector Ŝ at the head of each term.* + +### A13. Connected-cluster requirement +Every Hamiltonian fragment "must share at least one index with every cluster operator +on its right" (C&S p. 92); diagrams in which any T bar fails to contract with the +interaction are discarded. All T-vertex lines that do not reach the interaction must +be external. + +### A14. Kucharski–Bartlett sign sequences (enumeration bookkeeping) +Unique connectivities are enumerated by assigning "+" to particle and "−" to hole +lines below the Hamiltonian interaction line / above the cluster bars, and combining +the signs in all unique ways, e.g. the five unique sequences of (V_N T1 T1 T2)_c +in the T2 equation (C&S pp. 93–95; same idea as KST slides 17–18 "signatures": +T1 = +−, T2 = ++−−, V fragments = ++−−, etc.). + +--- + +## Part B — what CSV + DF adds (overlay annotations only) + +**B0. Skeleton invariance.** The CSV-CC diagram is a canonical Goldstone diagram +(Part A) plus annotations. Deleting all domain superscripts, C / CsC boxes and the +K label must leave a textbook diagram (e.g. R2 term 43 → the (V_N T1^4)_C doubles +diagram, C&S pp. 93–95). + +**B1. Domain superscripts.** A PNO virtual is written a^{ij} = virtual a restricted +to pair domain ⟨ij⟩ (SeQuant `a_1`; pao-ccsd-df.md Notation). Particle-line +labels always carry their domain. Occupied indices are ordinary MOs (no domain). + +**B2. DF split of the interaction (K).** Under density fitting +g_{pq,rs} = Σ_K B^K_{pq} B^K_{rs}: each half-vertex is one three-center factor +`g{bra;ket;Κ}` and **the dashed interaction line between the two half-vertices +carries the auxiliary label K** (write just "K"). The two `g` factors of a term are +identified by their shared Κ index. A term with `f{...}` instead has the one-body +×-vertex and no K. + +**B3. Slot → line mapping.** Every SeQuant tensor prints as `label{bra;ket[;aux]}`. +Apply A4 directly: **bra slot = outgoing line, ket slot = incoming line** at that +vertex. For amplitudes `t{virtuals;occupieds}`: virtuals are bra (created particles, +arrows leaving the bar upward), occupieds are ket (created holes: downward arrows +entering the bar from above). For the projector `Ŝ{occ;virt}` the occupieds are +bra-side (external hole lines flow from the top edge down into the diagram) and the +virtuals ket-side (external particle lines flow up out of the diagram). + +**B4. C annotation (PAO→PNO transform on an interaction leg).** g and f store their +*virtual* indices in the PAO basis (μ̃); a factor `C{μ̃;a}` (or transposed) +converts that leg to a PNO. Draw ONE fermion line for the chain g—C—a and put a +boxed **C** on the arrow near the g end: g —C→ a^{ij}. Occupied legs of g/f are MO +indices — never annotated. A leg with C is always a particle line. + +**B5. CsC annotation (domain-changing bridge).** A chain +`C{a;μ̃} · s{μ̃;μ̃'} · C{μ̃';a'}` (s = the PAO–PAO overlap metric) joins +two PNO indices of *different* domains into one particle line. Draw one line and put +a boxed **CsC** at mid-height. Typical use: connecting an external a_1^{i1 i2} to an +amplitude virtual a_5^{i3}. + +**B6. Fermion lines = connector chains.** Formally: C and s are not vertices. A +fermion line is a maximal chain of indices linked through C and s factors whose two +endpoints are slots of *terminal* tensors (g, f, t, Ŝ) — always exactly one bra slot +and one ket slot (this is the CSV image of A3/A4). Chain content ⇒ annotation: +one C → "C"; C·s·C → "CsC"; none → bare line. + +**B7. Amplitude vertex label.** `\hat{t}_{occ…}^{csv-virt…}` under the bar, e.g. +t{a_3;i_1} → t̂_{i1}^{a3^{i1}}; t{a_1,a_2;i_1,i_2} → +t̂_{i1 i2}^{a1^{i1 i2} a2^{i1 i2}}. Amplitudes are kept in CSV/PNO form — never +expanded into PAO (pao-ccsd-df.md Notation) — so C/CsC boxes never sit *between* a +T bar and its own index label; they sit on the line above it. + +**B8. External lines and Ŝ.** The external indices are exactly those listed in +`Ŝ{i…;a…}`; Ŝ is the biorthogonal particle-symmetrizer/projector marking +R = ⟨P|H̄|0⟩ — not a contracted tensor (pao-ccsd-df.md Notation). It replaces the +P(pq) machinery of A12 in the spin-traced biorthogonal formulation. + +**B9. Signs/weights.** The printed rational coefficient of each SeQuant term already +absorbs the (−1)^{h+l} sign (A9), the equivalent-line/vertex weights (A10/A11) and +the closed-shell spin-tracing factors; the diagram is drawn for structure, and the +A9–A11 count can be used as a consistency check (e.g. term 43: h = 4 hole lines, +l = 2 quasiloops pairing (i1,a1),(i2,a2) ⇒ (−1)^6 = +1, matching its coefficient +1). + +--- + +## Quick checklist for auditing a drawn CSV-CC diagram + +1. Every interaction half-vertex: exactly 1 in + 1 out (A3). +2. Every T1 bar: one particle up-arrow leaving + one hole down-arrow entering; + T_m bar: m of each (A2/A6). +3. Arrow directions: particles up, holes down, everywhere along a path (A1). +4. bra/ket of every SeQuant slot honored as out/in (A4/B3). +5. External lines = Ŝ indices only, run to the top edge, paired (i_k, a_k) + per the bra determinant (A7/B8). +6. Connectedness: every T bar shares ≥1 line with the interaction (A13). +7. K on the dashed line iff the term has two `g` factors sharing a Κ (B2). +8. C on every g/f particle leg; CsC on every C·s·C bridge; domains on every + particle label; no annotation on occupied lines (B1, B4, B5). +9. Strip annotations → a textbook Goldstone diagram (B0). diff --git a/utilities/diagram_gen/csv_diagram.py b/utilities/diagram_gen/csv_diagram.py new file mode 100644 index 0000000000..952cabebec --- /dev/null +++ b/utilities/diagram_gen/csv_diagram.py @@ -0,0 +1,1181 @@ +#!/usr/bin/env python3 +""" +csv_diagram.py — parse a SeQuant CSV-CCSD(-DF) term string (format of pao-ccsd-df.md) +and render the corresponding CSV-CC Goldstone diagram with matplotlib. + +Conventions implemented (see csv-cc-diagram-rules.md in this directory for citations): + * hole line = downward arrow (blue), occupied labels i_n + * particle line = upward arrow (red), PNO labels a_n^{dom} + * two-electron interaction = ONE vertex = two half-vertices (dots) joined by a + dashed line; under DF the dashed line carries the auxiliary label K + * one-body (Fock) interaction = a dot joined by a dashed line to a "×" marker + (A2: "the point of action ... marked by the interaction line —×") + * bra index = outgoing (creation); ket index = incoming (annihilation); + every half-vertex has exactly 1 incoming + 1 outgoing line + * fermion flow of a line: from the terminal node where its index sits in a BRA + slot to the terminal node where it sits in a KET slot + * amplitude vertex = solid horizontal bar at the bottom, label \\hat{t}_{occ}^{csv-virts} + * CSV overlay: PNO domain as superscript a^{ij}; PAO->PNO transform on a g/f leg + annotated "C"; a domain-changing bridge C·s·C annotated "CsC"; DF split of the + interaction annotated "K" on the dashed line + * external lines = indices listed in the projector Ŝ{occ;virt}; drawn to the top + +CSV OVERLAY KNOB (rule B0, skeleton invariance) + draw(..., show_csv=False) removes every CSV/DF annotation — domain superscripts, + C / CsC / s boxes, the K label and the g_L/g_R half-vertex tags — leaving the + plain canonical CC Goldstone diagram. show_aux= overrides the K/g_L/g_R + part alone (it defaults to show_csv). + +INPUT + Any SeQuant term string is accepted — inline, from a file, or on stdin. The + ':N-C-S' suffixes, markdown list numbering and backticks are stripped, so a + line copied straight out of pao-ccsd*.md works as-is: + + python3 csv_diagram.py '' out.png # inline + python3 csv_diagram.py --file term.txt out.png # from a file + sed -n '281p' ../pao-ccsdt-df.md | python3 csv_diagram.py - out.png + python3 csv_diagram.py --file many.txt outdir/ # N terms -> outdir/ + python3 render_all.py --src ../pao-ccsdtq-df.md --only R4 # whole section + +LAYOUT + * geometry is computed analytically (no renderer round-trip): text extents come + from a scratch canvas in DATA units, so the figure is built at exactly 1 data + unit = UNIT inches and every label is placed by a collision-avoiding greedy + search against vertices, bars, other labels and the fermion lines themselves. + * amplitude bars, half-vertices and external stubs are all centred on a common + axis, so mirror-image terms render as mirror images. + * two fermion lines joining the SAME pair of endpoints (e.g. the particle+hole + pair of a T1 hanging off a single interaction dot) are bowed apart into a + slim pointy ellipse (the canonical Goldstone bubble); lines that merely run + together (a shared apex, an external stub grazing a half-vertex) are fanned + apart by the same mechanism, so no two lines are ever drawn on top of one + another. + * C / CsC / s ride ON their fermion line and slide along it to dodge obstacles. + +VERIFIED ON + pao-ccsd-df.md (86), pao-ccsdt-df.md (201) and pao-ccsdtq-df.md (373 terms, + T1..T4 amplitude bars) — see audit_diagrams.py: no label/label, label/vertex or + line/line overlaps in any of them, in either variant. + +LIMITATIONS (not yet handled): + * line-crossing minimisation is only heuristic (amplitude/vertex ordering by + mean partner position); very crowded terms may still cross lines, which is + legal but not always pretty. Lines never coincide and labels never overlap, + but a label box may sit on top of one foreign line (1/86 CCSD, 8/201 CCSDT + and 32/373 CCSDTQ terms); it keeps an opaque background and stays readable. + * assumes at most ONE interaction (<=2 g half-vertices sharing one Κ, or one f). + True for every CCSD/CCSDT/CCSDTQ residual term in the files above. + * does not compute the diagram sign/weight (h, loops, equivalent lines); + it draws and annotates only. The printed coefficient is taken verbatim. + * an s{μ̃;μ̃} that is NOT flanked by two C's (never occurs in the files) would + be annotated as a bare "s". +""" + +import argparse +import os +import re +import sys +import unicodedata +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from matplotlib.lines import Line2D +import numpy as np + +HOLE = "#1f4fd8" # blue, arrow down +PART = "#c0281f" # red, arrow up +INT = "#555555" # gray dashed interaction + +Y_T, Y_HV, Y_EXT = 1.0, 3.9, 6.5 + +UNIT = 1.0 # inches per data unit => 72 pt per data unit +FS_IDX, FS_ANN, FS_T, FS_V, FS_TITLE, FS_LEG = 12, 11, 11, 12, 13, 10 + +APEX_STEP = 0.85 # spacing of the excitation apexes on one amplitude bar +BAR_HALF = 0.30 # bar overhang beyond the outermost apex +T_GAP = 0.75 # clear space between neighbouring amplitude bars/labels +HV_GAP = 3.4 # minimum separation of the two interaction half-vertices +EXT_BIAS = 0.14 # tie-break only: holes lean left of their anchor, particles right +EXT_GAP = 0.55 # clear space between neighbouring external labels +STUB = 1.15 # length of the Fock ---x stub + +LAST_AUDIT = None # geometry of the most recent draw() — used by audit_diagrams.py + +# ---------------------------------------------------------------------------- +# parsing +# ---------------------------------------------------------------------------- + +TENSOR_RE = re.compile(r"([^\s*(){}]+)\{([^{}]*)\}(?::[A-ZN\-]+)?") +COEFF_RE = re.compile(r"^\s*([+-]?\d+(?:/\d+)?)\s+") +MD_ITEM = re.compile(r"^\s*\d+\.\s+(.*)$") # "12. ``" markdown list line + + +def clean_term(s): + """Accept a raw SeQuant term as pasted from anywhere: a bare expression, a + markdown list item (``12. `...` ``), backtick-quoted, or with stray + surrounding whitespace / a leading '+'.""" + s = s.strip() + m = MD_ITEM.match(s) + if m and "{" in m.group(1): + s = m.group(1).strip() + s = s.strip("`").strip() + if s.startswith("+") and not s[1:2].isdigit(): + s = s[1:].strip() + return s + + +def _norm(s): + return unicodedata.normalize("NFC", s) + + +def kind_of(idx): + """occ / pno / pao / aux from the index name.""" + if idx.startswith("i"): + return "occ" + if idx.startswith("a"): + return "pno" + if "μ" in idx: + return "pao" + if "Κ" in idx or "K" in idx: + return "aux" + raise ValueError(f"unknown index kind: {idx!r}") + + +def split_outside(s, sep): + """split on sep, but not inside <...> (PNO domain brackets).""" + parts, depth, cur = [], 0, [] + for ch in s: + if ch == "<": + depth += 1 + elif ch == ">": + depth -= 1 + if ch == sep and depth == 0: + parts.append("".join(cur)); cur = [] + else: + cur.append(ch) + parts.append("".join(cur)) + return parts + + +def parse_index(tok): + """'a_3' -> ('a_3', ['i_1','i_2']); 'i_1' -> ('i_1', [])""" + tok = tok.strip() + m = re.match(r"([^<>]+)(?:<([^<>]*)>)?$", tok) + name = m.group(1).strip() + dom = [d.strip() for d in m.group(2).split(",")] if m.group(2) else [] + return name, dom + + +class Tensor: + def __init__(self, label, slots): + self.label = label # 'g','f','t','s','C','S^' + self.bra, self.ket = slots[0], slots[1] # lists of (name, dom) + self.aux = slots[2] if len(slots) > 2 else [] + + def __repr__(self): + return f"{self.label}{{{self.bra};{self.ket};{self.aux}}}" + + +def parse_term(term): + """Return (coeff:str, tensors:[Tensor]). Parens are irrelevant (associativity).""" + term = _norm(clean_term(term)) + if "{" not in term: + raise ValueError(f"no tensor found in input: {term[:80]!r}") + m = COEFF_RE.match(term) + coeff = m.group(1) if m else "1" + tensors = [] + for lab, body in TENSOR_RE.findall(term): + lab = _norm(lab) + if lab in ("Ŝ", "Ŝ", "S^"): + lab = "S^" + slots = [[parse_index(t) for t in split_outside(part, ",") if t.strip()] + for part in body.split(";")] + while len(slots) < 2: + slots.append([]) + tensors.append(Tensor(lab, slots)) + return coeff, tensors + + +# ---------------------------------------------------------------------------- +# graph construction: terminals + C/s connector chains -> fermion lines +# ---------------------------------------------------------------------------- + +TERMINALS = {"g", "f", "t", "S^"} + + +class Node: + """A drawable vertex: half-vertex ('hv'), amplitude ('t'), or external ('ext').""" + def __init__(self, kind, tag, tensor=None): + self.kind, self.tag, self.tensor = kind, tag, tensor + self.x = self.y = 0.0 + + def __repr__(self): + return f"<{self.kind}:{self.tag}>" + + +class FermionLine: + def __init__(self, src, dst, src_idx, dst_idx, nC, nS): + self.src, self.dst = src, dst # Node (bra end -> ket end) + self.src_idx, self.dst_idx = src_idx, dst_idx # (name, dom) at each end + self.nC, self.nS = nC, nS # of C / s tensors on the chain + # hole vs particle from terminal index kinds + kinds = {kind_of(src_idx[0]), kind_of(dst_idx[0])} + self.is_hole = "occ" in kinds + + @property + def annotation(self): + if self.nS: + return "CsC" if self.nC >= 2 else "sC" if self.nC else "s" + return "C" if self.nC == 1 else ("CC" if self.nC else None) + + +class UF: + def __init__(self): + self.p = {} + + def find(self, x): + self.p.setdefault(x, x) + while self.p[x] != x: + self.p[x] = self.p[self.p[x]] + x = self.p[x] + return x + + def union(self, a, b): + self.p[self.find(a)] = self.find(b) + + +def build_graph(tensors): + """Returns (nodes, lines, hvs, hv_pair_or_None, externals dict).""" + uf = UF() + conn_count = {} # root -> [nC, nS] (filled after unions) + connectors = [] # (label, p, q) + terminal_occ = {} # index name -> list of (node, 'bra'|'ket', (name,dom)) + nodes, hvs, tnodes = [], [], [] + S_hat = None + + for T in tensors: + if T.label in ("C", "s"): + p, q = T.bra[0][0], T.ket[0][0] + uf.union(p, q) + connectors.append((T.label, p, q)) + elif T.label in ("g", "f"): + n = Node("hv", f"{T.label}{len(hvs)}", T) + hvs.append(n); nodes.append(n) + elif T.label == "t": + n = Node("t", f"t{len(tnodes)}", T) + tnodes.append(n); nodes.append(n) + elif T.label == "S^": + S_hat = T + else: + raise ValueError(f"unhandled tensor {T}") + + def add_term(node, side, idx): + terminal_occ.setdefault(uf.find(idx[0]), []).append((node, side, idx)) + + for n in nodes: + for idx in n.tensor.bra: + add_term(n, "bra", idx) + for idx in n.tensor.ket: + add_term(n, "ket", idx) + + # externals: one ext node per Ŝ index. Ŝ bra=occ (creation at top edge), + # Ŝ ket=virt. Fermion flow: occ Ŝ->partner (down); virt partner->Ŝ (up). + ext_nodes = {} + if S_hat is not None: + for idx in S_hat.bra: # occupied externals + en = Node("ext", idx[0]); en.ext_idx = idx + ext_nodes[idx[0]] = en + add_term(en, "bra", idx) + for idx in S_hat.ket: # virtual externals + en = Node("ext", idx[0]); en.ext_idx = idx + ext_nodes[idx[0]] = en + add_term(en, "ket", idx) + + # count connectors per chain + for lab, p, q in connectors: + r = uf.find(p) + c = conn_count.setdefault(r, [0, 0]) + c[0 if lab == "C" else 1] += 1 + + lines = [] + for root, occs in terminal_occ.items(): + if len(occs) != 2: + raise ValueError(f"chain {root} has {len(occs)} terminals (want 2): {occs}") + (n1, s1, i1), (n2, s2, i2) = occs + if s1 == s2: + raise ValueError(f"chain {root}: both terminals are {s1} slots") + if s1 == "ket": + (n1, s1, i1), (n2, s2, i2) = (n2, s2, i2), (n1, s1, i1) + nC, nS = conn_count.get(root, [0, 0]) + lines.append(FermionLine(n1, n2, i1, i2, nC, nS)) + + # interaction pairing: g's sharing an aux index (DF), else lone f/g + hv_pair = None + if len(hvs) == 2: + hv_pair = (hvs[0], hvs[1]) + elif len(hvs) > 2: + raise ValueError(">2 interaction half-vertices not supported") + return nodes + list(ext_nodes.values()), lines, hvs, hv_pair, ext_nodes + + +# ---------------------------------------------------------------------------- +# labels +# ---------------------------------------------------------------------------- + +class Opts: + """Rendering knobs. show_csv=False => canonical CC Goldstone skeleton (rule B0). + + The three overlay layers can also be toggled independently: + csv C / CsC / s boxes on the fermion lines + dom PNO domain superscripts a^{ij} + aux the DF label K on the interaction line (+ the g_L/g_R half-vertex tags) + """ + def __init__(self, show_csv=True, show_aux=None, show_dom=None): + self.csv = bool(show_csv) + self.aux = self.csv if show_aux is None else bool(show_aux) + self.dom = self.csv if show_dom is None else bool(show_dom) + + +def _sub(name): + """'a_3' -> 'a_3' latex; 'μ̃_1' -> '\\tilde{\\mu}_1'""" + if "μ" in name: + n = name.split("_")[-1] + return rf"\tilde{{\mu}}_{{{n}}}" + base, _, n = name.partition("_") + return rf"{base}_{{{n}}}" if n else base + + +def idx_latex(idx, opts): + name, dom = idx + s = _sub(name) + if dom and opts.dom: + s += "^{\\," + "\\,".join(_sub(d) for d in dom) + "}" + return f"${s}$" + + +def t_label(tensor, opts): + occ = "\\,".join(_sub(n) for n, _ in tensor.ket) + vir = "\\,".join(_sub(n) + ("^{" + "\\,".join(_sub(d) for d in dom) + "}" + if (dom and opts.dom) else "") + for n, dom in tensor.bra) + return rf"$\hat{{t}}_{{\,{occ}}}^{{\,{vir}}}$" + + +# ---------------------------------------------------------------------------- +# text metrics (DATA units) — measured once on a scratch canvas +# ---------------------------------------------------------------------------- + +_SCRATCH = None +_SIZE_CACHE = {} + + +def text_size(s, fontsize): + """(width, height) of a rendered (math)text string, in data units.""" + key = (s, fontsize) + hit = _SIZE_CACHE.get(key) + if hit is not None: + return hit + global _SCRATCH + if _SCRATCH is None: + _SCRATCH = plt.figure(figsize=(2, 2), dpi=100) + _SCRATCH.canvas.draw() + art = _SCRATCH.text(0.5, 0.5, s, fontsize=fontsize) + bb = art.get_window_extent(_SCRATCH.canvas.get_renderer()) + art.remove() + res = (bb.width / _SCRATCH.dpi / UNIT, bb.height / _SCRATCH.dpi / UNIT) + _SIZE_CACHE[key] = res + return res + + +def box_of(x, y, w, h, ha="center", va="center", pad=0.06): + if ha == "center": + x0, x1 = x - w / 2, x + w / 2 + elif ha == "left": + x0, x1 = x, x + w + else: + x0, x1 = x - w, x + if va == "center": + y0, y1 = y - h / 2, y + h / 2 + elif va == "bottom": + y0, y1 = y, y + h + else: + y0, y1 = y - h, y + return (x0 - pad, y0 - pad, x1 + pad, y1 + pad) + + +# ---------------------------------------------------------------------------- +# geometry helpers +# ---------------------------------------------------------------------------- + +def _perp(u): + return np.array([u[1], -u[0]]) + + +def bow_path(p0, p1, bow, n=41): + """Quadratic-Bezier arc from p0 to p1 bulging by `bow` (perp, +=right of p0->p1). + Two opposite bows between the same endpoints make a pointy ellipse (lens).""" + p0 = np.asarray(p0, float); p1 = np.asarray(p1, float) + d = p1 - p0 + L = float(np.hypot(*d)) + if L < 1e-9 or abs(bow) < 1e-9: + return np.array([p0, p1]) + ctrl = (p0 + p1) / 2 + _perp(d / L) * 2.0 * bow + t = np.linspace(0, 1, n)[:, None] + return (1 - t) ** 2 * p0 + 2 * (1 - t) * t * ctrl + t ** 2 * p1 + + +def path_at(pts, f): + """point + unit tangent at arc-length fraction f of a polyline.""" + seg = np.diff(pts, axis=0) + ln = np.hypot(seg[:, 0], seg[:, 1]) + cum = np.concatenate([[0.0], np.cumsum(ln)]) + tot = cum[-1] + if tot < 1e-9: + return pts[0].copy(), np.array([0.0, 1.0]) + s = float(np.clip(f, 0.0, 1.0)) * tot + k = int(np.clip(np.searchsorted(cum, s) - 1, 0, len(seg) - 1)) + u = (s - cum[k]) / max(ln[k], 1e-12) + return pts[k] + seg[k] * u, seg[k] / max(ln[k], 1e-12) + + +def _overlap(a, b): + dx = min(a[2], b[2]) - max(a[0], b[0]) + dy = min(a[3], b[3]) - max(a[1], b[1]) + return dx * dy if (dx > 0 and dy > 0) else 0.0 + + +def _seg_hits_box(p, q, box): + """Liang-Barsky segment/AABB test.""" + x0, y0, x1, y1 = box + tmin, tmax = 0.0, 1.0 + dx, dy = q[0] - p[0], q[1] - p[1] + for pp, qq in ((-dx, p[0] - x0), (dx, x1 - p[0]), + (-dy, p[1] - y0), (dy, y1 - p[1])): + if abs(pp) < 1e-12: + if qq < 0: + return False + else: + r = qq / pp + if pp < 0: + if r > tmax: + return False + tmin = max(tmin, r) + else: + if r < tmin: + return False + tmax = min(tmax, r) + return tmin <= tmax + + +class Placer: + """Greedy collision-avoiding label placement against boxes and line segments.""" + + def __init__(self): + self.boxes = [] # occupied rectangles + self.segs = [] # (p, q, owner) + self.tags = [] # parallel to boxes: (tag, is_label) + + def block(self, box, tag="obstacle", is_label=False, owner=None): + self.boxes.append(box) + self.tags.append((tag, is_label, id(owner) if owner is not None else None)) + return box + + def block_text(self, x, y, s, fs, ha="center", va="center", pad=0.06): + w, h = text_size(s, fs) + return self.block(box_of(x, y, w, h, ha, va, pad)) + + def add_path(self, pts, owner, k=12): + pts = np.asarray(pts, float) + idx = np.unique(np.linspace(0, len(pts) - 1, min(k + 1, len(pts))).astype(int)) + for a, b in zip(idx[:-1], idx[1:]): + self.segs.append((pts[a], pts[b], owner)) + + def cost(self, box, owner=None): + """Any overlap with an existing box is disqualifying-expensive; sitting on + a foreign fermion line is merely ugly.""" + c = 0.0 + for b in self.boxes: + ov = _overlap(box, b) + if ov > 1e-9: + c += 3.0 + 25.0 * ov + for p, q, o in self.segs: + if _seg_hits_box(p, q, box): + c += 0.30 if o is owner else 0.85 + return c + + def place(self, cands, w, h, owner=None, pad=0.06, tag="label"): + """cands: [((x, y), preference_penalty)] — lower penalty = more wanted.""" + best = None + for x, y, pen in cands: + box = box_of(x, y, w, h, "center", "center", pad) + c = self.cost(box, owner) + pen + if best is None or c < best[0]: + best = (c, (x, y), box) + self.block(best[2], tag, True, owner) + return best[1] + + +# ---------------------------------------------------------------------------- +# layout +# ---------------------------------------------------------------------------- + +def _bar_half(tensor): + m = max(len(tensor.bra), len(tensor.ket), 1) + return APEX_STEP * (m - 1) / 2 + BAR_HALF + + +def assign_apexes(tnodes, lines): + """Map (t-node, line) -> apex point on the bar. The legs of one excitation + (particle+hole pair) MEET at a shared apex: a T_m vertex is m V-shapes.""" + slots = {} + for n in tnodes: + bra_names = [nm for nm, _ in n.tensor.bra] # virtuals -> particle legs + ket_names = [nm for nm, _ in n.tensor.ket] # occupieds -> hole legs + m = max(len(bra_names), len(ket_names), 1) + bra_line = [None] * len(bra_names) + ket_line = [None] * len(ket_names) + for L in lines: + if L.src is n: + idx = L.src_idx + elif L.dst is n: + idx = L.dst_idx + else: + continue + nm = idx[0] + if nm in bra_names: + bra_line[bra_names.index(nm)] = L + elif nm in ket_names: + ket_line[ket_names.index(nm)] = L + + def other_x(L): + if L is None: + return n.x + o = L.dst if L.src is n else L.src + return o.x + + pairs = [(bra_line[k] if k < len(bra_line) else None, + ket_line[k] if k < len(ket_line) else None) for k in range(m)] + pairs.sort(key=lambda pr: float(np.mean([other_x(L) for L in pr + if L is not None] or [n.x]))) + if m == 1: + apex_xs = [n.x] + else: + apex_xs = [n.x - APEX_STEP * (m - 1) / 2 + APEX_STEP * k for k in range(m)] + for k, (bl, kl) in enumerate(pairs): + for L in (bl, kl): + if L is not None: + slots[(n, id(L))] = (apex_xs[k], n.y) + for L in lines: # safety: stray leg -> centre + if (L.src is n or L.dst is n) and (n, id(L)) not in slots: + slots[(n, id(L))] = (n.x, n.y) + n.apex_xs = apex_xs + return slots + + +def _spread(pos, half, gap): + """Push overlapping items apart symmetrically (centroid preserving).""" + pos = list(pos) + for _ in range(400): + moved = False + for k in range(len(pos) - 1): + need = half[k] + half[k + 1] + gap + d = pos[k + 1] - pos[k] + if d < need - 1e-9: + push = (need - d) / 2.0 + pos[k] -= push + pos[k + 1] += push + moved = True + if not moved: + break + return pos + + +def _symmetrize(hvs): + if len(hvs) != 2: + return + a, b = hvs + if a.x > b.x: + a.x, b.x = b.x, a.x + mid = (a.x + b.x) / 2.0 + halfgap = max(HV_GAP, b.x - a.x) / 2.0 + a.x, b.x = mid - halfgap, mid + halfgap + + +def layout(nodes, lines, hvs, opts): + """Place every node; returns the apex map.""" + tnodes = [n for n in nodes if n.kind == "t"] + exts = [n for n in nodes if n.kind == "ext"] + + # which half-vertices does each amplitude touch? (left / straddling / right) + touch = {n: set() for n in tnodes} + for L in lines: + for a, b in ((L.src, L.dst), (L.dst, L.src)): + if a.kind == "t" and b.kind == "hv": + touch[a].add(b.tag) + + def group(n): + s = touch[n] + if len(hvs) == 2: + if s == {hvs[0].tag}: + return 0 + if s == {hvs[0].tag, hvs[1].tag}: + return 1 + if s == {hvs[1].tag}: + return 2 + return 1 + tnodes.sort(key=group) + + # --- bottom row: spacing driven by bar width AND amplitude-label width + half = [max(text_size(t_label(n.tensor, opts), FS_T)[0] / 2, _bar_half(n.tensor)) + for n in tnodes] + x = 0.0 + for k, n in enumerate(tnodes): + if k: + x += half[k - 1] + half[k] + T_GAP + n.x, n.y = x, Y_T + if tnodes: # centre the row on 0 + c = (tnodes[0].x + tnodes[-1].x) / 2 + for n in tnodes: + n.x -= c + + def hv_partner_xs(hv, slots=None): + out = [] + for L in lines: + for a, b in ((L.src, L.dst), (L.dst, L.src)): + if a is hv and b.kind == "t": + out.append(slots[(b, id(L))][0] if slots else b.x) + return out + + partner_of = {} + line_of = {} + for en in exts: + for L in lines: + if L.src is en or L.dst is en: + partner_of[en] = L.dst if L.src is en else L.src + line_of[en] = L + + def place_hvs(slots=None): + for j, hv in enumerate(hvs): + px = hv_partner_xs(hv, slots) + hv.x = float(np.mean(px)) if px else (HV_GAP * (j - 0.5)) + hv.y = Y_HV + _symmetrize(hvs) + + def ext_pref(en, slots=None): + p = partner_of.get(en) + if p is None: + return 0.0 + if p.kind == "t" and slots is not None: + px = slots[(p, id(line_of[en]))][0] + else: + px = p.x + return px + (-EXT_BIAS if kind_of(en.tag) == "occ" else EXT_BIAS) + + # provisional -> apexes -> refined half-vertices -> apexes -> externals + place_hvs() + for en in exts: + en.x, en.y = ext_pref(en), Y_EXT + slots = assign_apexes(tnodes, lines) + place_hvs(slots) + slots = assign_apexes(tnodes, lines) + + if exts: + exts.sort(key=lambda e: (ext_pref(e, slots), e.tag)) + halfw = [text_size(idx_latex(e.ext_idx, opts), FS_IDX)[0] / 2 for e in exts] + prefs = [ext_pref(e, slots) for e in exts] + for _ in range(3): + pos = _spread(prefs, halfw, EXT_GAP) + # keep external stubs from spearing an interaction dot + bumped = False + for k, en in enumerate(exts): + p = partner_of.get(en) + if p is None or p.kind == "hv": + continue + p0 = slots[(p, id(line_of[en]))] if p.kind == "t" else (p.x, p.y) + if abs(Y_EXT - p0[1]) < 1e-9: + continue + f = (Y_HV - p0[1]) / (Y_EXT - p0[1]) + xh = p0[0] + (pos[k] - p0[0]) * f + for hv in hvs: + if abs(xh - hv.x) < 0.45: + d = 0.45 - abs(xh - hv.x) + s = 1.0 if pos[k] >= hv.x else -1.0 + prefs[k] = pos[k] + s * (d / max(f, 0.2) + 0.15) + bumped = True + if not bumped: + break + pos = _spread(prefs, halfw, EXT_GAP) + for en, x in zip(exts, pos): + en.x, en.y = x, Y_EXT + + # --- global centring: mirror-image terms come out as mirror images + allx = [n.x for n in tnodes + exts + hvs] + for n in tnodes: + allx += list(getattr(n, "apex_xs", [])) + if allx: + c = (min(allx) + max(allx)) / 2.0 + for n in tnodes + exts + hvs: + n.x -= c + slots = assign_apexes(tnodes, lines) + return slots + + +# ---------------------------------------------------------------------------- +# drawing +# ---------------------------------------------------------------------------- + +def _sides(bow, p, q): + """Perpendicular sides to try first: the bulge side of a bowed line, else the + side pointing away from the diagram axis (keeps mirror terms mirror-imaged).""" + if abs(bow) > 1e-9: + return (1, -1) if bow > 0 else (-1, 1) + want = 1.0 if p[0] >= 0 else -1.0 + first = 1 if q[0] * want >= 0 else -1 + return (first, -first) + + +def _line_label_pairs(L, opts): + """(text_near_src, text_near_dst) — a bridge shows a different index at each + end, an ordinary line shows one label. PAO ends borrow the other end's name.""" + ks, kd = kind_of(L.src_idx[0]), kind_of(L.dst_idx[0]) + a = L.src_idx if ks != "pao" else (L.dst_idx if kd != "pao" else L.src_idx) + b = L.dst_idx if kd != "pao" else (L.src_idx if ks != "pao" else L.dst_idx) + ta, tb = idx_latex(a, opts), idx_latex(b, opts) + return (ta, tb) if ta != tb else (ta, None) + + +def draw(term, title=None, out="diagram.png", verbose=True, + show_csv=True, show_aux=None, show_dom=None, legend=True, dpi=170): + opts = Opts(show_csv, show_aux, show_dom) + coeff, tensors = parse_term(term) + nodes, lines, hvs, hv_pair, ext_nodes = build_graph(tensors) + slots = layout(nodes, lines, hvs, opts) + tnodes = [n for n in nodes if n.kind == "t"] + + if verbose: + print(f"coeff={coeff}; {len(hvs)} half-vertices; {len(tnodes)} amplitudes; " + f"{len(ext_nodes)} external lines") + for L in lines: + print(f" {'hole' if L.is_hole else 'part'} " + f"{L.src}({L.src_idx[0]}) -> {L.dst}({L.dst_idx[0]})" + f" ann={L.annotation}") + + def endpoint(node, L): + return np.array(slots[(node, id(L))] if node.kind == "t" else (node.x, node.y), + float) + + # ---- fermion-line geometry (bowing lines that share both endpoints) ----- + geo = {} + groups = {} + for L in lines: + p, q = endpoint(L.src, L), endpoint(L.dst, L) + key = tuple(sorted([tuple(np.round(p, 4)), tuple(np.round(q, 4))])) + groups.setdefault(key, []).append(L) + def _bow_line(L, bw): + """(re)build L's polyline as an arc bulging `bw` off the straight chord.""" + p, q = endpoint(L.src, L), endpoint(L.dst, L) + lo, hi = (p, q) if (p[1], p[0]) <= (q[1], q[0]) else (q, p) + pts = bow_path(lo, hi, bw) + if np.hypot(*(pts[0] - p)) > np.hypot(*(pts[-1] - p)): + pts = pts[::-1] + geo[id(L)] = dict(pts=pts, bow=bw, color=HOLE if L.is_hole else PART) + + # (a) lines sharing BOTH endpoints -> pointy ellipse (Goldstone bubble) + lens = set() + for key, grp in groups.items(): + n = len(grp) + if n == 1: + _bow_line(grp[0], 0.0) + continue + grp.sort(key=lambda L: (not L.is_hole, L.src_idx[0])) + chord = np.hypot(*(np.array(key[1]) - np.array(key[0]))) + span = float(np.clip(0.10 * chord, 0.14, 0.30)) # slim lens + for k, L in enumerate(grp): + _bow_line(L, -span + 2 * span * k / (n - 1)) + lens.add(id(L)) + + # (b) lines that merely RUN TOGETHER (sharing one endpoint, or an external + # stub grazing a half-vertex) are bowed apart too — just enough to read + # as two lines. Iterated, because bowing can uncover a new near-miss. + def _samples(pts, n=24): + d = np.cumsum(np.r_[0.0, np.hypot(*np.diff(pts, axis=0).T)]) + if d[-1] < 1e-9: + return np.repeat(pts[:1], n, axis=0) + u = np.linspace(0, d[-1], n) + return np.column_stack([np.interp(u, d, pts[:, 0]), + np.interp(u, d, pts[:, 1])]) + + def _runs_together(A, B, smp, hair=0.11, frac=0.30): + a, b = smp[id(A)], smp[id(B)] + d = np.hypot(a[:, None, 0] - b[None, :, 0], a[:, None, 1] - b[None, :, 1]) + return max((d.min(axis=1) < hair).mean(), (d.min(axis=0) < hair).mean()) > frac + + def _key(L): + return (not L.is_hole, L.src_idx[0], L.dst_idx[0]) + + # A bundle is fanned as a whole; the bundle graph only ever GROWS (so the + # fan never oscillates), and the fan widens only once membership is stable. + uf2, sep = UF(), 0.48 + for _ in range(6): + smp = {id(L): _samples(geo[id(L)]["pts"]) for L in lines} + hits = [(A, B) for i, A in enumerate(lines) for B in lines[i + 1:] + if _runs_together(A, B, smp)] + if not hits: + break + grew = False + for A, B in hits: + if uf2.find(id(A)) != uf2.find(id(B)): + uf2.union(id(A), id(B)) + grew = True + if not grew: + sep = min(sep * 1.3, 0.95) + comps = {} + for L in lines: + comps.setdefault(uf2.find(id(L)), []).append(L) + for grp in comps.values(): + if len(grp) < 2: + continue + free = [L for L in sorted(grp, key=_key) if id(L) not in lens] + if len(free) == 1: # push it off whatever it runs with + other = [L for L in grp if L is not free[0]] + b0 = np.sign(geo[id(other[0])]["bow"] or 1.0) if other else 1.0 + _bow_line(free[0], -b0 * sep) + elif free: + bows = [(k - (len(free) - 1) / 2) * sep for k in range(len(free))] + cap = max(abs(b) for b in bows) # keep arcs from wandering + if cap > 0.55: + bows = [b * 0.55 / cap for b in bows] + for L, bw in zip(free, bows): + _bow_line(L, bw) + + # ---- obstacles --------------------------------------------------------- + P = Placer() + texts = [] # (x, y, s, color, fs, ha, va, style) + + for n in tnodes: + x0, x1 = n.apex_xs[0] - BAR_HALF, n.apex_xs[-1] + BAR_HALF + P.block((x0 - 0.10, n.y - 0.14, x1 + 0.10, n.y + 0.14), "bar") + s = t_label(n.tensor, opts) + w, h = text_size(s, FS_T) + P.block(box_of(n.x, n.y - 0.34, w, h, "center", "top"), f"t:{s}", True) + texts.append((n.x, n.y - 0.34, s, "black", FS_T, "center", "top", "t")) + + for hv in hvs: + P.block((hv.x - 0.20, hv.y - 0.20, hv.x + 0.20, hv.y + 0.20), "dot") + + for en in [n for n in nodes if n.kind == "ext"]: + s = idx_latex(en.ext_idx, opts) + w, h = text_size(s, FS_IDX) + col = HOLE if kind_of(en.tag) == "occ" else PART + P.block(box_of(en.x, en.y + 0.22, w, h, "center", "bottom"), f"ext:{s}", True) + texts.append((en.x, en.y + 0.22, s, col, FS_IDX, "center", "bottom", "ext")) + + for L in lines: + P.add_path(geo[id(L)]["pts"], owner=L) + + # ---- interaction ------------------------------------------------------- + dash_seg = None + xmark = None + if len(hvs) == 2: + a, b = hvs + dash_seg = ((a.x, a.y), (b.x, b.y)) + P.segs.append((np.array(dash_seg[0]), np.array(dash_seg[1]), "dash")) + if opts.aux: + aux = a.tensor.aux[0][0] if a.tensor.aux else None + if aux: + w, h = text_size("$K$", FS_V) + # K belongs ON the dashed line; step sideways before stepping off it + cands = [((a.x + b.x) / 2 + d, a.y + dy, 0.30 * abs(d) + 1.30 * abs(dy)) + for d in (0, -0.55, 0.55, -1.1, 1.1, -1.65, 1.65) + for dy in (0, 0.42, -0.42)] + x, y = P.place(cands, w, h, owner="dash", tag="K") + texts.append((x, y, "$K$", INT, FS_V, "center", "center", "anno_int")) + for hv, lbl, side in ((a, "$g_L$", -1.0), (b, "$g_R$", 1.0)): + w, h = text_size(lbl, FS_V) + r = 0.30 + w / 2 + cands = [] + for k, (dx, dy) in enumerate(((side * r, 0.44), (side * r, -0.44), + (-side * r, 0.44), (-side * r, -0.44))): + for grow in (0.0, 0.45, 0.95): + cands.append((hv.x + dx + np.sign(dx) * grow, hv.y + dy, + 0.18 * k + 0.55 * grow)) + x, y = P.place(cands, w, h, owner="dash", tag=lbl) + texts.append((x, y, lbl, INT, FS_V, "center", "center", "hv")) + elif len(hvs) == 1: + hv = hvs[0] + best, bestc = (1.0, STUB), None + for s in (1.0, -1.0): + for ln in (STUB, STUB + 0.45, STUB + 0.9): + tip = hv.x + s * ln + box = (min(hv.x, tip) - 0.1, hv.y - 0.28, + max(hv.x, tip) + 0.24, hv.y + 0.28) + # the x marker itself must sit clear of everything + c = P.cost(box) + 3.0 * P.cost((tip - 0.26, hv.y - 0.26, + tip + 0.26, hv.y + 0.26)) + if bestc is None or c < bestc - 1e-9: + best, bestc = (s, ln), c + s, ln = best + dash_seg = ((hv.x, hv.y), (hv.x + s * ln, hv.y)) + xmark = (hv.x + s * ln, hv.y) + P.block((min(hv.x, xmark[0]) - 0.05, hv.y - 0.22, + max(hv.x, xmark[0]) + 0.22, hv.y + 0.22)) + + # ---- CSV annotations (C / CsC / s) ------------------------------------- + for L in lines: + ann = L.annotation if opts.csv else None + if not ann: + continue + g = geo[id(L)] + pts = g["pts"] + w, h = text_size(f"${ann}$", FS_ANN) + # "C" belongs next to the g/f end; bridges sit mid-line + # C / CsC / s ride ON the line (a bead threaded on the fermion line); + # they slide ALONG it to dodge obstacles, and only step aside as a + # last resort. "C" stays near its g/f end (rule B4). + if ann == "C" and (L.src.kind == "hv" or L.dst.kind == "hv"): + near_src = L.src.kind == "hv" + fr = [0.18, 0.26, 0.34, 0.42, 0.12, 0.50, 0.58, 0.66] + if not near_src: + fr = [1.0 - f for f in fr] + else: + fr = [0.50, 0.42, 0.58, 0.34, 0.66, 0.26, 0.74, 0.18, 0.82] + cands = [(*path_at(pts, f)[0], 0.12 * j) for j, f in enumerate(fr)] + for j, f in enumerate(fr[:4]): # last-resort step-aside + p, t = path_at(pts, f) + q = _perp(t) + off = abs(q[0]) * w / 2 + abs(q[1]) * h / 2 + 0.16 + for k, s in enumerate(_sides(g["bow"], p, q)): + cands.append((p[0] + s * q[0] * off, p[1] + s * q[1] * off, + 2.2 + 0.12 * j + 0.14 * k)) + x, y = P.place(cands, w, h, owner=L, tag=f"ann:{ann}") + texts.append((x, y, f"${ann}$", g["color"], FS_ANN, "center", "center", "anno")) + + # ---- index labels ------------------------------------------------------ + for L in lines: + g = geo[id(L)] + pts = g["pts"] + ta, tb = _line_label_pairs(L, opts) + jobs = [] + src_ext = L.src.kind == "ext" + dst_ext = L.dst.kind == "ext" + if tb is None: # one name for the line + if not (src_ext or dst_ext): # externals label at the top + jobs.append((ta, [0.38, 0.52, 0.26, 0.64, 0.14, 0.76, 0.88])) + else: # bridge: name at each end + if not src_ext: + jobs.append((ta, [0.24, 0.34, 0.16, 0.44, 0.08, 0.52])) + if not dst_ext: + jobs.append((tb, [0.76, 0.66, 0.84, 0.56, 0.92, 0.48])) + for s, fr in jobs: + w, h = text_size(s, FS_IDX) + cands = [] + for extra in (0.17, 0.52, 0.95, 1.55): # progressively roomier + for j, f in enumerate(fr): + p, t = path_at(pts, f) + q = _perp(t) + off = abs(q[0]) * w / 2 + abs(q[1]) * h / 2 + extra + for k, sg in enumerate(_sides(g["bow"], p, q)): + cands.append((p[0] + sg * q[0] * off, p[1] + sg * q[1] * off, + 1.40 * (extra - 0.17) + 0.09 * j + 0.13 * k)) + x, y = P.place(cands, w, h, owner=L, tag=f"idx:{s}") + texts.append((x, y, s, g["color"], FS_IDX, "center", "center", "idx")) + + # ---- audit record (see audit_diagrams.py) ------------------------------ + global LAST_AUDIT + LAST_AUDIT = dict(boxes=list(P.boxes), tags=list(P.tags), + segs=[(tuple(p), tuple(q), o) for p, q, o in P.segs]) + + # ---- extents ----------------------------------------------------------- + xs, ys = [], [] + for L in lines: + pts = geo[id(L)]["pts"] + xs += [pts[:, 0].min(), pts[:, 0].max()] + ys += [pts[:, 1].min(), pts[:, 1].max()] + for b in P.boxes: + xs += [b[0], b[2]]; ys += [b[1], b[3]] + if dash_seg: + xs += [dash_seg[0][0], dash_seg[1][0]]; ys += [dash_seg[0][1], dash_seg[1][1]] + xs = xs or [-1.0, 1.0]; ys = ys or [0.0, 1.0] + m = 0.30 + xr = max(abs(min(xs) - m), abs(max(xs) + m)) + xmin, xmax = -xr, xr # symmetric frame about the diagram axis + ymin, ymax = min(ys) - m, max(ys) + m + if ext_nodes: + ymax = max(ymax, Y_EXT + 0.60) + + # ---- figure ------------------------------------------------------------ + LM = RM = 0.14 + TM = 0.52 if title else 0.14 + BM = 0.62 if legend else 0.14 + figw = (xmax - xmin) * UNIT + LM + RM + figh = (ymax - ymin) * UNIT + TM + BM + + need = 0.6 + if title: + need = max(need, text_size(title, FS_TITLE)[0] * UNIT + 0.4) + if legend: + need = max(need, 6.1) + if figw < need: # widen symmetrically + extra = (need - figw) / UNIT + xmin -= extra / 2; xmax += extra / 2 + figw = (xmax - xmin) * UNIT + LM + RM + + fig = plt.figure(figsize=(figw, figh), dpi=dpi) + ax = fig.add_axes([LM / figw, BM / figh, + (xmax - xmin) * UNIT / figw, (ymax - ymin) * UNIT / figh]) + ax.set_xlim(xmin, xmax); ax.set_ylim(ymin, ymax) + ax.axis("off") + + if ext_nodes: # projection boundary (no caption) + ax.axhline(Y_EXT, color="black", lw=0.6, ls=(0, (1, 4)), alpha=0.35, zorder=1) + + # amplitude bars + for n in tnodes: + ax.plot([n.apex_xs[0] - BAR_HALF, n.apex_xs[-1] + BAR_HALF], [n.y, n.y], + color="black", lw=4, solid_capstyle="butt", zorder=5) + + # interaction + if dash_seg: + ax.plot([dash_seg[0][0], dash_seg[1][0]], [dash_seg[0][1], dash_seg[1][1]], + color=INT, lw=1.6, ls=(0, (6, 4)), zorder=2) + for hv in hvs: + ax.plot([hv.x], [hv.y], "o", color=INT, ms=11, zorder=6) + if xmark is not None: # one-body vertex: ---x + ax.plot([xmark[0]], [xmark[1]], marker="x", color=INT, ms=13, mew=2.6, + zorder=6, ls="none") + + # fermion lines + arrows + for L in lines: + g = geo[id(L)] + pts = g["pts"] + ax.plot(pts[:, 0], pts[:, 1], color=g["color"], lw=1.8, zorder=3, + solid_capstyle="round") + p, t = path_at(pts, 0.5) + ax.annotate("", xy=p + t * 0.28, xytext=p - t * 0.28, + arrowprops=dict(arrowstyle="-|>", color=g["color"], lw=0, + mutation_scale=22), zorder=6) + + # labels + for x, y, s, col, fs, ha, va, style in texts: + if style in ("anno", "anno_int"): + ax.text(x, y, s, ha=ha, va=va, fontsize=fs, color=col, zorder=9, + bbox=dict(boxstyle="round,pad=0.16", fc="#fff7d6", ec=col, + alpha=0.96)) + elif style == "t": + ax.text(x, y, s, ha=ha, va=va, fontsize=fs, color=col, + fontweight="bold", zorder=9) + elif style == "hv": + ax.text(x, y, s, ha=ha, va=va, fontsize=fs, color=col, + fontweight="bold", zorder=9) + else: + ax.text(x, y, s, ha=ha, va=va, fontsize=fs, color=col, zorder=8, + bbox=dict(boxstyle="round,pad=0.10", fc="white", ec="none", + alpha=0.80)) + + if legend: + if len(hvs) == 2: + ilab = "two-body vertex" + (" ($K$)" if opts.aux else "") + ihandle = Line2D([0], [0], color=INT, lw=1.6, ls="--") + else: + ilab = "one-body vertex ($-\\!-\\!\\times$)" + ihandle = Line2D([0], [0], color=INT, lw=1.6, ls="--") + leg = [Line2D([0], [0], color=PART, lw=2, label="particle line (up)"), + Line2D([0], [0], color=HOLE, lw=2, label="hole line (down)"), + (ihandle, ilab)] + fig.legend(handles=[leg[0], leg[1], leg[2][0]], + labels=[leg[0].get_label(), leg[1].get_label(), leg[2][1]], + loc="lower center", ncol=3, frameon=False, fontsize=FS_LEG, + bbox_to_anchor=(0.5, 0.004)) + if title: + ax.set_title(f"{title} (coefficient {coeff})", fontsize=FS_TITLE, pad=8) + + out = os.path.expanduser(out) + fig.savefig(out, dpi=dpi) + plt.close(fig) + if verbose: + print("saved:", out) + return out + + +# ---------------------------------------------------------------------------- +# demo: R2 terms from pao-ccsd-df.md +# ---------------------------------------------------------------------------- + +TERM43 = ("Ŝ{i_1,i_2;a_1,a_2}:N-C-S * g{i_3;μ̃_1;Κ_1}:N-C-S * " + "C{μ̃_1;a_3}:N-C-S * g{i_4;μ̃_2;Κ_1}:N-C-S * C{μ̃_2;a_4}:N-C-S * " + "s{μ̃_3;μ̃_4}:N-C-S * C{a_1;μ̃_3}:N-C-S * C{μ̃_4;a_5}:N-C-S * " + "s{μ̃_5;μ̃_6}:N-C-S * C{a_2;μ̃_5}:N-C-S * C{μ̃_6;a_6}:N-C-S * " + "t{a_3;i_1}:N-C-S * t{a_4;i_2}:N-C-S * t{a_5;i_3}:N-C-S * " + "t{a_6;i_4}:N-C-S") + +# R2 term 30: (oo|oo) interaction — both half-vertices carry two HOLE legs; +# both external particles reach the amplitudes only through CsC bridges. +TERM30 = ("((Ŝ{i_1,i_2;a_1,a_2}:N-C-S * ((s{μ̃_1;μ̃_2}:N-C-S * " + "(C{μ̃_2;a_3}:N-C-S * t{a_3;i_3}:N-C-S)) * " + "C{a_1;μ̃_1}:N-C-S)) * (g{i_3;i_1;Κ_1}:N-C-S * " + "g{i_4;i_2;Κ_1}:N-C-S)) * ((s{μ̃_3;μ̃_4}:N-C-S * (C{μ̃_4;a_4}:N-C-S * " + "t{a_4;i_4}:N-C-S)) * C{a_2;μ̃_3}:N-C-S)") + +# R2 term 49: (vv|vv) interaction + one T2 — each half-vertex has particle-in +# (from the T2, via C) and particle-out (external, via C); the T2 also carries +# both external hole lines directly. +TERM49 = ("(Ŝ{i_1,i_2;a_1,a_2}:N-C-S * (((g{μ̃_1;μ̃_2;Κ_1}:N-C-S * " + "C{a_1;μ̃_1}:N-C-S) * C{μ̃_2;a_3}:N-C-S) * " + "((g{μ̃_3;μ̃_4;Κ_1}:N-C-S * C{a_2;μ̃_3}:N-C-S) * " + "C{μ̃_4;a_4}:N-C-S))) * t{a_3,a_4;i_1,i_2}:N-C-S") + + +def _cli(): + ap = argparse.ArgumentParser( + description="render CSV-CC Goldstone diagram(s) from SeQuant term string(s)", + epilog="the term may be given inline, with --file, or on stdin; markdown " + "list numbering and backticks are stripped automatically") + ap.add_argument("term", nargs="?", help="SeQuant term string ('-' = read stdin)") + ap.add_argument("out", nargs="?", default="diagram.png", + help="output PNG (or output DIRECTORY when --file holds many terms)") + ap.add_argument("--file", metavar="PATH", + help="read the term(s) from a file ('-' = stdin); every line " + "holding a tensor expression is rendered") + ap.add_argument("--title", default=None) + ap.add_argument("--no-csv", action="store_true", + help="canonical CC skeleton: drop C/CsC/s, domains and K") + ap.add_argument("--no-dom", action="store_true", help="drop domain superscripts only") + ap.add_argument("--no-aux", action="store_true", help="drop the K / g_L,g_R tags only") + ap.add_argument("--no-legend", action="store_true") + args = ap.parse_args() + kw = dict(show_csv=not args.no_csv, + show_dom=False if args.no_dom else None, + show_aux=False if args.no_aux else None, + legend=not args.no_legend) + + src = args.file or ("-" if args.term == "-" else None) + if src is not None: # file / stdin: 1..N terms + # with --file the single positional is the OUTPUT ("--file in.txt out.png") + out = args.out if (args.out != "diagram.png" or not args.file) else \ + (args.term or args.out) + text = sys.stdin.read() if src == "-" else open(src, encoding="utf-8").read() + terms = [ln for ln in text.splitlines() if "{" in ln and "}" in ln] + if not terms: + ap.error(f"no tensor expression found in {src}") + if len(terms) == 1: + draw(terms[0], title=args.title, out=out, **kw) + else: + outdir = out if (out.endswith("/") or "." not in os.path.basename(out)) \ + else (os.path.dirname(out) or ".") + os.makedirs(outdir, exist_ok=True) + print(f"{len(terms)} terms -> {outdir}/") + for k, t in enumerate(terms, 1): + draw(t, title=(args.title or "term") + f" {k}", verbose=False, + out=os.path.join(outdir, f"term{k:03d}.png"), **kw) + elif args.term: + draw(args.term, title=args.title, out=args.out, **kw) + else: + print("no term given — drawing the three built-in demo terms here:") + for nm, tm, ti in (("term43", TERM43, r"CSV-CCSD $R_2$ term 43"), + ("term30", TERM30, r"CSV-CCSD $R_2$ term 30 (oo|oo)"), + ("term49", TERM49, r"CSV-CCSD $R_2$ term 49 (vv|vv)")): + draw(tm, title=ti, out=f"csv_ccsd_{nm}.png", verbose=False, **kw) + + +if __name__ == "__main__": + _cli() diff --git a/utilities/diagram_gen/examples/ccsd_E3_fock_bubble.png b/utilities/diagram_gen/examples/ccsd_E3_fock_bubble.png new file mode 100644 index 0000000000..8c0115b811 Binary files /dev/null and b/utilities/diagram_gen/examples/ccsd_E3_fock_bubble.png differ diff --git a/utilities/diagram_gen/examples/ccsd_R2t30_oooo_csc_bridge.png b/utilities/diagram_gen/examples/ccsd_R2t30_oooo_csc_bridge.png new file mode 100644 index 0000000000..f4c988ee02 Binary files /dev/null and b/utilities/diagram_gen/examples/ccsd_R2t30_oooo_csc_bridge.png differ diff --git a/utilities/diagram_gen/examples/ccsd_R2t49_canonical.png b/utilities/diagram_gen/examples/ccsd_R2t49_canonical.png new file mode 100644 index 0000000000..f584c20f1f Binary files /dev/null and b/utilities/diagram_gen/examples/ccsd_R2t49_canonical.png differ diff --git a/utilities/diagram_gen/examples/ccsd_R2t49_vvvv_ladder.png b/utilities/diagram_gen/examples/ccsd_R2t49_vvvv_ladder.png new file mode 100644 index 0000000000..6df6eddaea Binary files /dev/null and b/utilities/diagram_gen/examples/ccsd_R2t49_vvvv_ladder.png differ diff --git a/utilities/diagram_gen/examples/ccsdt_R1t1_with_T3.png b/utilities/diagram_gen/examples/ccsdt_R1t1_with_T3.png new file mode 100644 index 0000000000..b1b94094ad Binary files /dev/null and b/utilities/diagram_gen/examples/ccsdt_R1t1_with_T3.png differ diff --git a/utilities/diagram_gen/render_all.py b/utilities/diagram_gen/render_all.py new file mode 100644 index 0000000000..8a767a6fb9 --- /dev/null +++ b/utilities/diagram_gen/render_all.py @@ -0,0 +1,167 @@ +#!/usr/bin/env python3 +""" +render_all.py — batch-render every term of a SeQuant equation file (sections +E / R1 / R2 / R3 / R4 of pao-ccsd-df.md, pao-ccsdt-df.md, pao-ccsdtq-df.md, ...) +to a CSV-CC Goldstone diagram PNG using the csv_diagram renderer. + +Two variants (rule B0, skeleton invariance): + --variant csv (default) full CSV/DF overlay: domains, C/CsC boxes, K + --variant canonical plain CC Goldstone skeleton (no C/s, no K) + --variant both render both sets + + python3 render_all.py # CCSD, CSV overlay + python3 render_all.py --variant both + python3 render_all.py --src ../pao-ccsdt-df.md --only R3 + python3 render_all.py --src ../pao-ccsdtq-df.md --only R4:1,2,3 --outdir /tmp/q + +Outputs go to diagrams//
_term.png (zero-padded); the default +outdir is derived from the source file name and the variant. A manifest INDEX.md +is written summarizing successes/failures per section. +""" + +import argparse +import os +import re + +from csv_diagram import draw # parse_term is used internally by draw + +HERE = os.path.dirname(os.path.abspath(__file__)) + +# Default equation file: the CCSD one, wherever this script happens to live. +SRC_CANDIDATES = [ + os.path.join(HERE, "pao-ccsd-df.md"), # alongside the script + os.path.join(os.path.dirname(HERE), "pao-ccsd-df.md"), # one directory up +] +SRC = next((p for p in SRC_CANDIDATES if os.path.exists(p)), SRC_CANDIDATES[0]) + +# "## E — ...", "## R1 — ...", ... "## R4 — ..." (CCSD / CCSDT / CCSDTQ files) +SECTION_RE = re.compile(r"^##\s+(E|R\d+)\b") +TERM_RE = re.compile(r"^\s*(\d+)\.\s+`(.+?)`\s*$") + + +def extract_terms(path): + """Return {section_tag: [(num, expr), ...]} in file order.""" + with open(path, encoding="utf-8") as fh: + lines = fh.readlines() + + sections = {} + current = None + for line in lines: + stripped = line.rstrip("\n") + if stripped.startswith("## "): + m = SECTION_RE.match(stripped) + current = m.group(1) if m else None # Notation etc. -> no section + if current: + sections.setdefault(current, []) + continue + if current is None: + continue + m = TERM_RE.match(stripped) + if m: + sections[current].append((int(m.group(1)), m.group(2))) + return sections + + +def section_order(sections): + """E first, then R1, R2, R3, ... in numeric order.""" + return ([t for t in ("E",) if t in sections] + + sorted((t for t in sections if t != "E"), key=lambda t: int(t[1:]))) + + +def render_set(sections, outdir, show_csv, kind="CSV-CCSD"): + os.makedirs(outdir, exist_ok=True) + results = {} # tag -> list of (num, filename_or_None, error_or_None) + for tag in section_order(sections): + terms = sections.get(tag, []) + results[tag] = [] + for num, expr in terms: + fname = f"{tag}_term{num:02d}.png" + out = os.path.join(outdir, fname) + title = f"{kind} {tag} term {num}" + try: + draw(expr, title=title, out=out, verbose=False, show_csv=show_csv) + results[tag].append((num, fname, None)) + except Exception as exc: # noqa: BLE001 - want to keep going + msg = f"{type(exc).__name__}: {exc}" + results[tag].append((num, None, msg)) + print(f"FAILED {tag} term {num}: {msg}") + write_index(results, outdir, show_csv, kind) + return results + + +def write_index(results, outdir, show_csv, kind="CSV-CCSD"): + head = f"# {kind} diagram batch — INDEX" + if not show_csv: + head += " (canonical skeleton)" + lines = [head, ""] + grand_ok = grand_fail = 0 + for tag in section_order(results): + rows = results.get(tag, []) + ok = [r for r in rows if r[2] is None] + fail = [r for r in rows if r[2] is not None] + grand_ok += len(ok) + grand_fail += len(fail) + lines.append(f"## {tag} ({len(ok)} rendered / {len(fail)} failed of {len(rows)})") + lines.append("") + for num, fname, err in rows: + if err is None: + lines.append(f"- term {num}: `{fname}`") + else: + lines.append(f"- term {num}: FAILED: {err}") + lines.append("") + lines.insert(1, "") + lines.insert(2, f"**Totals:** {grand_ok} rendered / {grand_fail} failed " + f"(of {grand_ok + grand_fail} terms).") + with open(os.path.join(outdir, "INDEX.md"), "w", encoding="utf-8") as fh: + fh.write("\n".join(lines) + "\n") + print(f"{outdir}: {grand_ok} rendered / {grand_fail} failed") + + +def filter_sections(sections, only): + """only = 'R2' or 'R2:16' or 'R2:16,49'""" + tag, _, nums = only.partition(":") + keep = {int(n) for n in nums.split(",") if n} if nums else None + return {t: [(n, e) for n, e in v if keep is None or n in keep] + for t, v in sections.items() if t == tag} + + +def main(): + ap = argparse.ArgumentParser(description=__doc__.splitlines()[1]) + ap.add_argument("--src", default=SRC, help="SeQuant equation markdown file") + ap.add_argument("--variant", choices=("csv", "canonical", "both"), default="csv") + ap.add_argument("--outdir", default=None, + help="output directory (single-variant runs only)") + ap.add_argument("--only", default=None, + help="restrict to e.g. R2 or R2:16 or R2:16,49") + ap.add_argument("--kind", default=None, + help="title prefix (default: derived from --src)") + args = ap.parse_args() + + src = os.path.abspath(os.path.expanduser(args.src)) + if not os.path.exists(src): + ap.error("equation file not found: %s\nlooked for a default in:\n %s\n" + "pass one explicitly with --src " + % (src, "\n ".join(SRC_CANDIDATES))) + sections = extract_terms(src) + if args.only: + sections = filter_sections(sections, args.only) + + stem = os.path.splitext(os.path.basename(src))[0] # e.g. pao-ccsdt-df + theory = stem.replace("pao-", "").replace("-df", "").upper() # CCSDT + kind = args.kind or f"CSV-{theory}" + base = "ccsd_all" if stem == "pao-ccsd-df" else f"{theory.lower()}_all" + + out = {} + if args.variant in ("csv", "both"): + d = (args.outdir if (args.outdir and args.variant == "csv") + else os.path.join(HERE, base)) + out["csv"] = render_set(sections, d, True, kind) + if args.variant in ("canonical", "both"): + d = (args.outdir if (args.outdir and args.variant == "canonical") + else os.path.join(HERE, base + "_canonical")) + out["canonical"] = render_set(sections, d, False, theory) + return out + + +if __name__ == "__main__": + main() diff --git a/utilities/diagram_gen/sample_terms.txt b/utilities/diagram_gen/sample_terms.txt new file mode 100644 index 0000000000..0541ea0150 --- /dev/null +++ b/utilities/diagram_gen/sample_terms.txt @@ -0,0 +1,10 @@ +# One SeQuant term per line; '#' lines are ignored (no braces). +# Rendered with: python3 csv_diagram.py --file sample_terms.txt out/ +# CCSD E term 3 — T1 + Fock vertex (bubble) +2 (f{i_1;μ̃_1}:N-C-S * C{μ̃_1;a_1}:N-C-S) * t{a_1;i_1}:N-C-S +# CCSD R2 term 49 — (vv|vv) ladder + T2 +(Ŝ{i_1,i_2;a_1,a_2}:N-C-S * (((g{μ̃_1;μ̃_2;Κ_1}:N-C-S * C{a_1;μ̃_1}:N-C-S) * C{μ̃_2;a_3}:N-C-S) * ((g{μ̃_3;μ̃_4;Κ_1}:N-C-S * C{a_2;μ̃_3}:N-C-S) * C{μ̃_4;a_4}:N-C-S))) * t{a_3,a_4;i_1,i_2}:N-C-S +# CCSD R2 term 30 — (oo|oo) + two T1 via CsC bridge +((Ŝ{i_1,i_2;a_1,a_2}:N-C-S * ((s{μ̃_1;μ̃_2}:N-C-S * (C{μ̃_2;a_3}:N-C-S * t{a_3;i_3}:N-C-S)) * C{a_1;μ̃_1}:N-C-S)) * (g{i_3;i_1;Κ_1}:N-C-S * g{i_4;i_2;Κ_1}:N-C-S)) * ((s{μ̃_3;μ̃_4}:N-C-S * (C{μ̃_4;a_4}:N-C-S * t{a_4;i_4}:N-C-S)) * C{a_2;μ̃_3}:N-C-S) +# CCSDT R1 term 1 — T3 amplitude +Ŝ{i_1;a_1}:N-C-S * ((((((g{i_2;μ̃_1;Κ_1}:N-C-S * C{μ̃_1;a_2}:N-C-S) * (g{i_3;μ̃_2;Κ_1}:N-C-S * C{μ̃_2;a_3}:N-C-S)) * t{a_2,a_3,a_4;i_3,i_1,i_2}:N-C-S) * C{μ̃_3;a_4}:N-C-S) * s{μ̃_4;μ̃_3}:N-C-S) * C{a_1;μ̃_4}:N-C-S)