diff --git a/.gitignore b/.gitignore index e601b6422..dc672ea0b 100644 --- a/.gitignore +++ b/.gitignore @@ -79,3 +79,6 @@ nsqso_born.inc docs/build/ tests/input_files/IOTestsComparison_BackUp/ tests/input_files/IOTestsComparison/**/*.BackUp + +# output of ./tests/test_manager.py +UNITTEST_proc/ diff --git a/aloha/template_files/aloha_functions.f b/aloha/template_files/aloha_functions.f index 01aa9d8ad..09dbd86a3 100644 --- a/aloha/template_files/aloha_functions.f +++ b/aloha/template_files/aloha_functions.f @@ -1923,3 +1923,43 @@ subroutine CombineAmpS(nb, ihels, iwfcts, W1, Wall, Amp) enddo return end + + subroutine sumw_1(w1, w2, wout) +c +c Sum two currents standing for the same off shell line: the four +c gluon current and the pair of three gluon vertices it factorises +c into carry the same colour factor, so the amplitude reading the sum +c gets both contributions from a single call. See +c HelasMatrixElement.get_quartic_current_sums. +c +c The two share their momentum, so only the wavefunction is added and +c everything else is taken over from the first one. +c + use ALOHA_OBJECT + implicit none + type(aloha) w1 + type(aloha) w2 + type(aloha) wout + + wout = w1 + wout%W(:) = w1%W(:) + w2%W(:) + + return + end + + + subroutine subw_1(w1, w2, wout) +c +c As sumw_1, for the contributions which enter with a minus sign. +c + use ALOHA_OBJECT + implicit none + type(aloha) w1 + type(aloha) w2 + type(aloha) wout + + wout = w1 + wout%W(:) = w1%W(:) - w2%W(:) + + return + end diff --git a/aloha/template_files/aloha_functions_fd.f b/aloha/template_files/aloha_functions_fd.f index 2696c80d8..fcac1ac46 100644 --- a/aloha/template_files/aloha_functions_fd.f +++ b/aloha/template_files/aloha_functions_fd.f @@ -2240,3 +2240,43 @@ subroutine define_gauge_dir(q, n) end + + subroutine sumw_1(w1, w2, wout) +c +c Sum two currents standing for the same off shell line: the four +c gluon current and the pair of three gluon vertices it factorises +c into carry the same colour factor, so the amplitude reading the sum +c gets both contributions from a single call. See +c HelasMatrixElement.get_quartic_current_sums. +c +c The two share their momentum, so only the wavefunction is added and +c everything else is taken over from the first one. +c + use ALOHA_OBJECT + implicit none + type(aloha) w1 + type(aloha) w2 + type(aloha) wout + + wout = w1 + wout%W(:) = w1%W(:) + w2%W(:) + + return + end + + + subroutine subw_1(w1, w2, wout) +c +c As sumw_1, for the contributions which enter with a minus sign. +c + use ALOHA_OBJECT + implicit none + type(aloha) w1 + type(aloha) w2 + type(aloha) wout + + wout = w1 + wout%W(:) = w1%W(:) - w2%W(:) + + return + end diff --git a/aloha/template_files/aloha_functions_loop.f b/aloha/template_files/aloha_functions_loop.f index 46561ce7d..8d558bbe0 100644 --- a/aloha/template_files/aloha_functions_loop.f +++ b/aloha/template_files/aloha_functions_loop.f @@ -3042,3 +3042,43 @@ subroutine olxxxx(p,ffmass,nhel,nsf,fo) c return end + + subroutine sumw_1(w1, w2, wout) +c +c Sum two currents standing for the same off shell line: the four +c gluon current and the pair of three gluon vertices it factorises +c into carry the same colour factor, so the amplitude reading the sum +c gets both contributions from a single call. See +c HelasMatrixElement.get_quartic_current_sums. +c +c The two share their momentum, so only the wavefunction is added and +c everything else is taken over from the first one. +c + use ALOHA_OBJECT + implicit none + type(aloha) w1 + type(aloha) w2 + type(aloha) wout + + wout = w1 + wout%W(:) = w1%W(:) + w2%W(:) + + return + end + + + subroutine subw_1(w1, w2, wout) +c +c As sumw_1, for the contributions which enter with a minus sign. +c + use ALOHA_OBJECT + implicit none + type(aloha) w1 + type(aloha) w2 + type(aloha) wout + + wout = w1 + wout%W(:) = w1%W(:) - w2%W(:) + + return + end diff --git a/docs/gluon-quartic-plan.md b/docs/gluon-quartic-plan.md new file mode 100644 index 000000000..2a8cb6d1c --- /dev/null +++ b/docs/gluon-quartic-plan.md @@ -0,0 +1,1020 @@ +# Pure-gluon amplitude optimisation — plan + +Branch `claude/gluon-amplitude-optimization-8706f5`. Everything is behind +`set merge_quartic_vertices` (off by default), so nothing changes until it is +set. It takes four values: + +| | | +|---|---| +| `False` | off, the default | +| `speed` | the current sums, and the diagram order which allows them | +| `slots` | the order which keeps fewest currents alive; no current sums | +| `auto` | `slots` when the matrix elements go to a gpu backend, `speed` otherwise, decided per output | + +It has to be a `set` option and not an `output` one, because the diagram order +is fixed while the diagrams are generated -- by `output` time it is already +too late. The interface pushes it onto `madgraph.merge_quartic_vertices`, +which is what the generation and the exporters read. + +`auto` is the exception, and only because `slots` is the `speed` order +reversed: it generates the `speed` order and +`MadGraphCmd.apply_quartic_diagram_order` reverses it at `output` time, once +the backend about to be handed the matrix elements is known. The choice is +made from the *matrix element* exporter rather than the output format, so +`output madevent --me_exporter=` -- one output feeding two +backends off a single `_curr_matrix_elements` -- follows the gpu. + +**speed against slots.** The wavefunction store is a stack frame on cpu and is +per thread on gpu, where at seven gluons it is about 24 kB a thread and caps +occupancy. So the trade goes opposite ways: + +| `g g > 5 g` | amplitude calls | slots | per thread | +|---|---|---|---| +| off | 7245 | 268 | 25.1 kB | +| `speed` | 6813 | 259 | 24.3 kB | +| `slots` | 7245 | **199** | **18.7 kB** | + +6% more arithmetic for 23% less memory. Measured on cpu `speed` wins; which +one a gpu wants has *not* been measured -- there is no device here. + +## Goal + +Each colour structure of the 4-gluon vertex carries the *same colour factor* +as the diagram obtained by splitting that vertex into two cubic ones. So the +quartic current and the cubic current carrying that colour factor can be +summed into one current, and the whole subtree below it emitted once instead +of twice — a Berends-Giele style recursion, without colour ordering. + +``` +TMP = VVVVk_1(W_a,W_b,W_c) + VVV1P0_1(VVV1P0_1(W_a,W_b), W_c) +``` +Both carry the same `1/P^2` through the same propagator, so this is a plain +sum. Same for amplitudes. + +## Established facts (measured, not assumed) + +**1. The colour identity holds exactly.** Grouping every amplitude piece by +its colour vector over the ColorBasis: + +| process | amplitude pieces | colour groups | cubic diagrams per group | +|---|---|---|---| +| `g g > g g` | 6 | 3 | exactly 1 | +| `g g > g g g` | 45 | 15 | exactly 1 | +| `g g > g g g g` | 510 | 105 | exactly 1 | + +Every quartic piece is colour-proportional (±1) to exactly one cubic diagram. +Merged count == pure-cubic diagram count == `(2n-5)!!`. `g g > 4g` also +verified in generated Fortran: folding 405 of 510 amplitudes and zeroing the +sources leaves `|M|^2` bit-identical. + +**2. The seed rule — forbid two 3-gluon vertices from sharing a line.** + +| process | full | seed | seed % | reconstructed by unrolling | missing | +|---|---|---|---|---|---| +| `g g > g g` | 4 | 1 | 25% | 4 | 0 | +| `g g > g g g` | 25 | 10 | 40% | 25 | 0 | +| `g g > g g g g` | 220 | 55 | 25% | 220 | 0 | +| `g g > g g g g g` | 2485 | 385 | 15.5% | 2485 | 0 | + +This is forced, not tuned. Unrolling a quartic vertex always yields two +**adjacent** cubic vertices (joined by the line that replaced the vertex), so +a diagram is reachable iff it has an adjacent cubic pair to contract back. +The diagrams with no such pair are exactly the ones that must be in the seed, +and every diagram either has such a pair or is in the seed — necessary and +sufficient, hence exact coverage. + +Note "no 3-gluon vertices **at all**" is the over-restrictive special case: it +loses the diagrams whose cubic vertices are non-adjacent (60 of 220 at six +gluons). + +**3. Reconstruction gives matched rootings for free.** A reconstructed diagram +and its quartic partner come from the *same* seed diagram, so they share a +decomposition and every current except the one being summed. This is the +property the whole optimisation needs. + +## Already committed + +| commit | what | +|---|---| +| `fcd8218b6` | link map (`unroll_quartic_vertices`, verified vs colour algebra, N=2..5) + amplitude sums | +| `5c9cdb301` | why the current sums fail on the un-rewired DAG | +| `1c1722ae8` | revert of the auxiliary-particle generation | +| `98d288f41` | `reroot_diagram`, validated 755/755 — probably NOT needed under the seed rule | +| `25fc6d1cc` | step 1, the seed rule inside `reduce_leglist` | +| `7da06bac7` | step 2+3, `expand_seed_diagrams` and the recorded links | +| `8e634cf9a` | step 4, `TMP = W1 + c*W4` at the amplitudes | + +Useful pieces to keep: `get_unrollable_quartic_vertices`, +`unroll_quartic_vertices`, `diagram_colour_signature`, `UnrollDiagramTag`, +`split_quartic_vertex`, `unrolled_diagram`, `get_quartic_amplitude_merges`, +`get_amplitude_merge_lines`. + +## Plan + +**Step 1 — enforce the seed rule in generation.** DONE, `25fc6d1cc`. Rejected +inside `reduce_leglist`, tracking by leg number which lines a cubic vertex +produced. Measured against the full generation filtered by an independently +written adjacency detector, comparing the diagrams and not only the counts: +1 / 10 / 55 / 385, and 4165 of 34300 at eight gluons. + +The closing vertex needs two cases, not one. A real n->0 interaction takes +the legs left over as lines coming in, but a 2->2 like reduction ends on the +*identity* vertex, which only states that its two legs are the two ends of +one line — and that line joins the two vertices which produced them. Missing +that left `g g > g g` with all four of its diagrams. + +**Step 2 — reconstruct the full set by unrolling the seed.** DONE, +`1b9474f69`, `expand_seed_diagrams`. Diagram count exactly the baseline +(4/25/220/2485), same diagrams by tag, `|M|^2` bit-identical at four and five +gluons and to 1e-14 at six and seven. + +The dedup has to run on the *glued* form. While the identity vertex is still +there the same diagram has several spellings, and a quartic vertex sitting +just in front of it is not yet the last one — which is what decides how ALOHA +indexes its colour structures. Without gluing first, the dedup caught nothing +(40 diagrams for `g g > g g g`) and 8 of 30 recorded colour chains disagreed +with the colour algebra. + +**Step 3 — record the link during reconstruction.** DONE, same commit, +`get_quartic_unroll_links`. 3 / 30 / 405 / 6300 links, every one agreeing +target for target with `unroll_quartic_vertices`, which stays as the +independent colour-algebra cross-check. + +**Step 4 — the current sum.** DONE at the amplitudes, `8e634cf9a`, and +blocked deeper. `TMP = W1 + c*W4` where a quartic current and its cubic +partner feed the same vertex; the amplitude reads the sum and the quartic +amplitude is never computed: + +``` +W(20) = W(19) +W(20)%W(:) = W(19)%W(:) + W(12)%W(:) +CALL VVV1_0(W(4),W(5),W(20),GC_10,AMP(33)) +``` + +The sum is shared by every amplitude reading it, so it pays for as many calls +as it has users: 60 amplitude calls for 30 sums at six gluons, 432 for 60 at +seven. Substituting several mothers of one amplitude also produces the +amplitude with all of them substituted, so every subset has to be a merge into +that same target weighing the product of the coefficients — that check is what +keeps the count honest. + +Two things it costs. `reuse_outdated_wavefunctions` works out when a slot is +free from the diagrams alone and cannot see the extra read, so it has to be +told, or the two currents get handed the same slot (`W(11) + W(11)`). And the +sums take a slot each, never reused: `NWAVEFUNCS` 51 -> 91 at six gluons. + +| | `g g > g g g` | `g g > g g g g` | `g g > 5 g` | +|---|---|---|---| +| helas calls | 94 -> 93 + 7 sums | 637 -> 612 + 30 | 8159 -> 7784 + 60 | +| JAMP lines | 131 -> 101 | 1082 -> 688 | 23672 -> 7864 | +| per call | 34.88 -> 35.12 s | 47.77 -> 46.02 s | 42.16 -> 39.80 s | + +**The sum stays at the amplitudes.** Measured +on the reconstructed matrix element (`g g > g g g g`): of the 275 places where +a cubic current is fed by another cubic current *and* the quartic partner +taking the same four lines exists, 225 are the last vertex — the amplitude +sum, which pitfall 7 says buys nothing — and 50 are genuine currents. **None +of the 50 pass the 1:1 consumer test.** At seven gluons, none of 135. Only at +five gluons do all 7 pass. + +Why, from a failing pair: quartic current 30 has consumers +`{Wav(3,30), Amp(3,4,6,30)x3, Amp(6,8,30), Amp(4,10,30)}` while its cubic +partner 53 has `{Amp(6,8,53), Amp(4,10,53), Amp(3,12,53)}`. Two correspond; +the rest do not, because the same diagram is rooted differently on the two +sides. + +That is forced, not a bug in the reconstruction. Expanding the seed produces +more (seed, choice) instances than there are diagrams, so some diagrams are +reached from several seeds: + +| process | seeds | (seed, choice) instances | diagrams | +|---|---|---|---| +| `g g > g g` | 1 | 4 | 4 | +| `g g > g g g` | 10 | 40 | 25 | +| `g g > g g g g` | 55 | 340 | 220 | +| `g g > 5 g` | 385 | 4900 | 2485 | + +A diagram carries one decomposition, so it can be rooted to match at most one +of its quartic partners, while the current sum needs the match at *every* +node. Fact 3 above holds per (seed, choice) and breaks under the dedup that +fact 2 requires. Not fixable by choosing the spelling more cleverly: the +counting alone rules it out, and `g g > g g` — the one row with no collision +— is also the one process where the seed rule reaches every partner. + +That is why the sum is taken only where the vertex reading it is an +*amplitude*: there is nothing above it, so no consumer can be handed a term +it must not have and the substitution is a plain swap of one argument. A sum +at a current would have to be carried through every current above it, and +those are shared. What would work there is a partial rewrite — hand the sum +only to the consumers which do correspond and leave W1 and W4 serving the +rest — which splits shared consumers and cascades upward. That is a DAG +rewriting problem, not this plan. + +The same counting is what leaves the two-substitution cases at seven gluons +on the table: their target has a source with both mothers substituted, but +the source with only the *other* one substituted is spelled with a different +rooting, so the subset check refuses it. 432 of the 864 amplitude calls the +structure allows — and the rest are not simply waiting to be picked up, see +"where to go next". + +**Step 5 — validate and time.** DONE. `|M|^2` for `g g > N g`, N=2..5, and +per-call timing from the shipped `check` driver, which already loops +`SMATRIX` when given a second argument (`./check 1000 20000`). + +| | `g g > g g g g` | `g g > 5 g` | +|---|---|---| +| flag off | 47.73 / 47.77 / 47.80 s | 42.15 / 42.16 s | +| flag on, before steps 1-4 | 47.76 s | 40.86 s | +| flag on, steps 1-3 only | 48.03 / 48.07 s | 40.41 s | +| flag on, at HEAD | 46.14 / 45.90 s | 39.80 s | + +and the code that produces it: + +| | helas calls | JAMP lines | +|---|---|---| +| flag off | 637 / 8159 | 1082 / 23672 | +| flag on, before steps 1-4 | 637 / 8159 | 688 / 8012 | +| flag on, steps 1-3 only | 672 / 8216 | 688 / 7864 | +| flag on, at HEAD | 612 + 30 / 7784 + 60 | 688 / 7864 | + +So the flag is worth **+3.8% at six gluons and +5.6% at seven**. The JAMP +shrink from `fcd8218b6` carries seven gluons on its own; six gluons only +turns positive with the current sum, because the reconstruction deviates from +the canonical decomposition — which is the whole point — and thereby weakens +the wavefunction CSE, 35 calls at six gluons and 57 at seven. Five gluons is +a wash (34.88 -> 35.12 s): 7 sums against 7 calls is too little to pay for +the 14 extra slots. + +`|M|^2` bit-identical at four and five gluons, 1e-15 at six and seven, and +unchanged for `g g > t t~ g g` and `u u~ > g g g`. + +With the flag off, `matrix.f` is byte-identical to `3b3ed9e85` for N=2..5. + +## Step 6 — madevent, and what it does to AMP2 + +`29b0c670e`. Two things had to give before the optimisation survived the +madevent path, neither of them about AMP2: + +1. **The sum has to be a `CALL`.** Helicity recycling rebuilds the whole DAG + from the calls alone, so a bare assignment was invisible to it — the summed + current never entered the graph and the amplitude reading it died on a + `KeyError`. Written as `CALL SUMW_1(W(a),W(b),W(c))` it is an ordinary + internal wavefunction with two mothers and everything downstream handles + it. `sumw_1`/`subw_1` live in `aloha_functions.f`; the coefficient is + restricted to ±1, which is all it has ever been. +2. **`hel_recycle.add_indices` could not index a statement-initial `AMP(`.** + The pattern ate the character in front of it and there is none at the + start of a line, so `AMP(31) = AMP(31) + AMP(1)` came out as + `AMP(31) = AMP( K,31) + AMP( K,1)`. Latent until now: nothing had ever + emitted a line beginning with `AMP(`. + +**AMP2 is left alone and picks up the merged amplitude**, which is the right +thing and turns out to matter more than the matrix element speedup. The fold +lines run before the AMP2 block, so the channel weight is +`|AMP_cubic + AMP_quartic|^2`. Nothing else was needed: `get_amp2_lines` +already skips any diagram with a four point vertex, so the folded amplitudes +were never referenced. + +That skip is the point. With the flag off, at six gluons: + +| | amplitudes computed | distinct AMP reaching AMP2 | +|---|---|---| +| flag off | 510 | **105** | +| flag on | 450 (+30 sums) | 105, now carrying all 510 | + +so four fifths of the amplitude — every quartic contribution — used to enter +**no** channel weight at all. Folding puts each of them in the channel whose +colour factor it shares, which is exactly where it belongs. + +**The integration is not measurably better or worse.** That is the answer, and +it took some care to get to, because `generate_events` wall time is a bad +metric here: the refine stage adapts, and the *same* configuration swings by a +factor four between seeds (217 s to 877 s with the flag off). Every number +below is mean +- standard error over independent seeds. + +`g g > g g g`, `generate_events`, 10000 events, three seeds: + +| | cross section | rel. error | ME cpu | +|---|---|---|---| +| flag off | 3.680-3.694e+07 pb | 0.327% | 52.7 s | +| flag on | 3.684-3.694e+07 pb | 0.300% | 49.0 s | + +`g g > g g g g`, `generate_events`, 2000 events, four seeds, and a survey-only +run — the same fixed number of points both ways, so the error measures the +channel weights and nothing else — over six seeds: + +| | rel. error (survey) | rel. error (full) | cpu (full) | +|---|---|---|---| +| flag off | 2.14% | 0.498% | 583 s | +| flag on | 1.93% | 0.473% | 580 s | + +Every difference is favourable and none of them is significant: + +| | off - on | | +|---|---|---| +| `g g > g g g` error | -8% | 1.8 sigma | +| `g g > g g g` cpu | -7% | 1.1 sigma | +| `g g > g g g g` survey error | -10% | 1.1 sigma | +| `g g > g g g g` full error | -5% | 0.8 sigma | +| `g g > g g g g` full cpu | -1% | 0.0 sigma | + +Cross sections agree everywhere. So the honest reading is that handing the +quartic contributions to the channel that shares their colour factor does not +hurt the integration, and may help it slightly — but the seed to seed scatter +is far larger than the effect, and nothing here is established at more than +two sigma. A campaign of tens of seeds would be needed to call it either way. + +## Step 7 — madmatrix + +`e93dfa1b1`. `SUMW_1`/`SUBW_1` as C++ templates next to `ALOHAOBJ` in +`cpp_hel_amps_h.inc` (they cannot use the `INLINE` macro, which the ALOHA +generated block defines further down the header), and the madmatrix writer +emits them exactly as the Fortran one does. + +**The colour amplitudes had to be sorted out first, and that was a live bug.** +`get_color_amplitudes` dropped every merge source from the JAMPs on the +assumption that the caller writes the amplitude sums to put them back. Only +the Fortran writer does, so C++ and python output with the flag set +was quietly losing four fifths of the amplitude. It now takes +`merge_quartic_amplitudes`; a backend which writes no sums keeps those +amplitudes in the JAMPs, where their own colour coefficients give the +identical result. So madmatrix gets the current sums, which really do remove +work, and leaves the amplitude level merges alone — there is no `AMP` array +there to fold into anyway, each amplitude going straight into the JAMPs. + +One bug in the shared writer surfaced here: a wavefunction number can be +listed by more than one diagram in the madmatrix matrix element (two objects +for the same current with the mothers ordered differently), so a sum was +written twice — 50 lines for 30 sums at six gluons. Harmless numerically, both +write the same value to the same slot, but wasted. Both writers now write each +sum once. + +|M|^2 to 1e-14 at five and six gluons (`FPTYPE=d`; the default mixed +precision build rounds the two to the same value), `CPPProcess.cc` +byte-identical with the flag off. + +**Speed is mixed, and worse than Fortran** (after step 8): + +| | amp calls | nwf | evt/s (sse4, FPTYPE=d) | | +|---|---|---|---|---| +| `g g > g g g` | 45 -> 38 | 12 -> 19 | 72820 -> 68010 | **-6.6%** | +| `g g > g g g g` | 510 -> 450 | 51 -> 86 | 2696 -> 2748 | **+1.9%** | + +There is no JAMP fold here to pay for the extra wavefunctions, so five gluons +loses outright: 7 sums against 7 saved amplitude calls does not cover a +wavefunction array half again as large. Note the slot count with the flag on +is worse in madmatrix than in Fortran (86 against 66 at six gluons), because +its matrix element carries duplicate wavefunctions -- two objects for the same +current with the mothers ordered differently, which MG5 does not merge. + +## Step 8 — recycle the slots the sums take + +`5b421bc80`. A sum used to get a slot of its own at the end, never reused, +which more than doubled `NWAVEFUNCS`. But a sum is an ordinary producer -- +written as soon as the later of its two currents is made, dead after the last +amplitude reading it -- so it goes through `reuse_outdated_wavefunctions` with +everything else. That also lets the cubic current die at the sum rather than +at the amplitude, since the amplitude no longer reads it. + +| | flag off | own slots | recycled | +|---|---|---|---| +| `g g > g g g` | 12 | 26 | **19** | +| `g g > g g g g` | 51 | 91 | **66** | +| `g g > g g g g` (madmatrix) | 51 | 111 | **86** | + +**A bug in `reuse_outdated_wavefunctions` had to be fixed first, and it was +not introduced here.** The same wavefunction can be listed by more than one +diagram in the madmatrix matrix element, and the allocator handed it a slot +again on the second listing, leaking the first. Harmless while the sums had +their own slots; once they shared the pool, `g g > g g g g` came out 0.2% +wrong. A wavefunction now takes one slot at its first appearance and keeps it +until its last use. Nothing moves with the flag off -- `matrix.f` and +`CPPProcess.cc` are byte-identical for N=2..4 in both backends. + +Worth it for madmatrix (-8.4% -> -6.6% at five gluons, +0.9% -> +1.9% at six) +and neutral for Fortran (+3.8% -> +3.9%): there the slot count was never the +bottleneck. The madevent run is unchanged, same cross section and error. + +## Results + +Everything below is `g g > N g` with the flag off against `speed`, on the +same machine. Standalone Fortran is the shipped `check` driver looping +`SMATRIX`; madmatrix is `check_sa.exe perf` built `FPTYPE=d` on `cppsse4` +(the default mixed precision build rounds the two to the same value and would +hide any difference). Two runs each, reproducible to about 0.1%. + +**Speed** + +| | standalone | | madmatrix | | +|---|---|---|---|---| +| | per call | | evt/s | | +| `g g > g g` | 11.00 -> 11.04 s | -0.4% | 875150 -> 878724 | +0.4% | +| `g g > g g g` | 34.90 -> 34.99 s | -0.3% | 72359 -> 66757 | **-7.7%** | +| `g g > g g g g` | 47.35 -> 45.61 s | **+3.7%** | 2699 -> 2784 | **+3.1%** | +| `g g > 5 g` | 43.05 -> 39.98 s | **+7.1%** | 41.75 -> 43.33 | **+3.8%** | + +Four gluons is a wash on both (there is nothing to sum: the only quartic +vertex is the whole amplitude). Five gluons loses on madmatrix, where the +wavefunction store grows by half and there is no JAMP fold to pay for it. Six +and seven gluons win on both, and the gain grows with the multiplicity. + +**Where it came from.** Three states: the flag off (what an unoptimised build +gives), the branch as it stood before this work (`3b3ed9e85`, only the +amplitude merges of `fcd8218b6`), and the flag on today. + +standalone Fortran, `slots / wavefunction calls + sums / amplitude calls`: + +| | flag off | before this work | at PR open | now | +|---|---|---|---|---| +| `g g > g g` | 5 / 7 / 6 | 5 / 7 / 6 | 5 / 7 / 6 | 5 / 7 / 6 | +| `g g > g g g` | 12 / 33 / 45 | 12 / 33 / 45 | 19 / 39+7 / 38 | 19 / 39+7 / 38 | +| `g g > g g g g` | 51 / 111 / 510 | 51 / 111 / 510 | 66 / 146+30 / 450 | **78 / 126+30 / 450** | +| `g g > 5 g` | 268 / 898 / 7245 | 268 / 898 / 7245 | 290 / 955+60 / 6813 | **259 / 925+60 / 6813** | + +madmatrix, `slots / amplitude calls`: + +| | flag off | before this work | at PR open | now | +|---|---|---|---|---| +| `g g > g g` | 5 / 6 | *wrong* | 5 / 6 | 5 / 6 | +| `g g > g g g` | 12 / 45 | *wrong* | 19 / 38 | 19 / 38 | +| `g g > g g g g` | 51 / 510 | *wrong* | 86 / 450 | **78 / 450** | +| `g g > 5 g` | 268 / 7245 | *wrong* | 320 / 6813 | **259 / 6813** | + +"wrong" is not a figure of speech: at `3b3ed9e85` `get_color_amplitudes` dropped +every merge source from the JAMPs unconditionally, and only the Fortran writer +wrote the sums putting them back, so with the flag set any C++ or python output +lost 405 of its 510 amplitudes at six gluons, silently. That is fixed here. + +Per-call time, standalone Fortran (lower is better) and madmatrix in evt/s +(higher is better): + +| | off | before | at PR open | now | +|---|---|---|---|---| +| fortran `g g > g g g g` | 47.77 s | 48.12 s | 45.61 s | **45.25 s (+5.3%)** | +| fortran `g g > 5 g` | 42.19 s | 40.88 s | 39.98 s | **40.10 s (+4.9%)** | +| madmatrix `g g > g g g` | 74229 | *wrong* | 66757 | 67566 (-9.0%) | +| madmatrix `g g > g g g g` | 2651 | *wrong* | 2784 | **2880 (+8.6%)** | + +Timings taken in batches, and there is about 1% of drift between batches, so +each column should be read against the `off` measured with it. + +**Memory — the wavefunction store**, which is what the optimisation costs. +`NWAVEFUNCS` in Fortran, `nwf` in madmatrix; bytes are 100 per wavefunction in +Fortran (4 complex, 4 reals, one int) and 192 in madmatrix on sse4 in double +(4 complex plus a 4-momentum, over a 2 event vector). + +| | off | on, slot each | on, recycled | | +|---|---|---|---|---| +| `g g > g g` | 5 | 5 | 5 | 500 B | +| `g g > g g g` | 12 | 26 | **19** | 1900 B (+58%) | +| `g g > g g g g` | 51 | 91 | **66** | 6600 B (+29%) | +| `g g > 5 g` | 268 | 321 | **290** | 29000 B (+8%) | + +madmatrix carries duplicate wavefunctions of its own, so its count with the +flag on is higher: 19 / 86 / 320 at five, six and seven gluons, against +19 / 66 / 290 in Fortran. The relative cost falls as the multiplicity rises +(+58% at five gluons, +19% at seven), which is why the trade turns positive +from six gluons on. + +**Work done per call** + +| | standalone helas calls | JAMP lines | madmatrix amplitude calls | jamp lines | +|---|---|---|---|---| +| `g g > g g` | 29 -> 29 | 23 -> 20 | 6 -> 6 | 34 -> 34 | +| `g g > g g g` | 94 -> 100 (+7 sums) | 131 -> 101 | 45 -> 38 | 370 -> 314 | +| `g g > g g g g` | 637 -> 642 (+30) | 1082 -> 688 | 510 -> 450 | 8170 -> 7210 | +| `g g > 5 g` | 8159 -> 7844 (+60) | 23672 -> 7864 | 7245 -> 6813 | 231850 -> 218026 | + +`|M|^2` is bit-identical at four and five gluons and agrees to 1e-14 at six +and seven, in both backends. With the flag off, `matrix.f` and `CPPProcess.cc` +are byte-identical to before any of this. + +## Full sweep -- speed, memory and generation time + +All three modes against the flag off, on the same machine, for both series. +Timings are the minimum of five runs of the shipped `check` driver looping +`SMATRIX` on a fixed phase space point; the minimum matters, because a single +run carries about 3% of noise and most of the effects here are smaller than +that. The floor of the method is 0.6%, measured on `g g > t t~`, where `off` +and `speed` produce a byte-identical `matrix.f` and still time 4.000 against +3.988 us. `|M|^2` agrees to 8.5e-15 or better on every row. + +*slots* is `NWAVEFUNCS`, the length of the wavefunction array, `TYPE(ALOHA) +W(NWAVEFUNCS)` -- not the number of wavefunctions computed, since +`reuse_outdated_wavefunctions` frees an entry as soon as its last reader has +run (898 wavefunction calls live in 268 slots at seven gluons). One entry is +104 bytes, measured with `storage_size`: four `complex*16`, `P(0:3)` and +`flv_index`, padded. One amplitude is 16 bytes. + +**`g g > N g`** + +| process | mode | generate | matrix.f | matrix.o | W slots | W array | AMP entries | amps computed | AMP array | per call | speed | +|---|---|---|---|---|---|---|---|---|---|---|---| +| g g > 2g | off | 1.8 s | 27 kB | 16 kB | 5 | 0.5 kB | 6 | 6 | 0.1 kB | 5.48 us | - | +| | speed | 1.5 s | 27 kB | 17 kB | 5 | 0.5 kB | 4 | 6 | 0.1 kB | 5.52 us | -1% | +| | slots | 1.6 s | 27 kB | 17 kB | 5 | 0.5 kB | 4 | 6 | 0.1 kB | 5.55 us | -1% | +| g g > 3g | off | 2.1 s | 40 kB | 28 kB | 12 | 1.2 kB | 45 | 45 | 0.7 kB | 87.50 us | - | +| | speed | 1.8 s | 39 kB | 29 kB | 19 | 1.9 kB | 24 | 38 | 0.4 kB | 90.25 us | -3% | +| | slots | 2.4 s | 40 kB | 29 kB | 12 | 1.2 kB | 16 | 45 | 0.2 kB | 93.00 us | -6% | +| g g > 4g | off | 2.8 s | 181 kB | 141 kB | 51 | 5.2 kB | 510 | 510 | 8.0 kB | 2.37 ms | - | +| | speed | 2.4 s | 165 kB | 146 kB | 78 | 7.9 kB | 316 | 450 | 4.9 kB | 2.29 ms | +4% | +| | slots | 2.4 s | 168 kB | 152 kB | 54 | 5.5 kB | 106 | 510 | 1.7 kB | 2.41 ms | -2% | +| g g > 5g | off | 35.2 s | 3.3 MB | 5.8 MB | 268 | 27.2 kB | 7245 | 7245 | 113.2 kB | 141.00 ms | - | +| | speed | 16.6 s | 2.5 MB | 3.4 MB | 259 | 26.3 kB | 5869 | 6813 | 91.7 kB | 132.75 ms | +6% | +| | slots | 18.4 s | 2.5 MB | 3.4 MB | 199 | 20.2 kB | 946 | 7245 | 14.8 kB | 134.25 ms | +5% | + +**`g g > t t~ N g`** + +| process | mode | generate | matrix.f | matrix.o | W slots | W array | AMP entries | amps computed | AMP array | per call | speed | +|---|---|---|---|---|---|---|---|---|---|---|---| +| g g > t t~ | off | 1.5 s | 26 kB | 16 kB | 5 | 0.5 kB | 3 | 3 | 0.0 kB | 4.00 us | - | +| | speed | 1.5 s | 26 kB | 16 kB | 5 | 0.5 kB | 3 | 3 | 0.0 kB | 3.99 us | +0% | +| | slots | 2.0 s | 26 kB | 16 kB | 5 | 0.5 kB | 3 | 3 | 0.0 kB | 4.09 us | -2% | +| g g > t t~ g | off | 1.8 s | 31 kB | 20 kB | 12 | 1.2 kB | 18 | 18 | 0.3 kB | 30.00 us | - | +| | speed | 1.6 s | 31 kB | 20 kB | 12 | 1.2 kB | 15 | 15 | 0.2 kB | 29.58 us | +1% | +| | slots | 1.7 s | 31 kB | 21 kB | 12 | 1.2 kB | 15 | 18 | 0.2 kB | 30.33 us | -1% | +| g g > t t~ 2g | off | 2.6 s | 68 kB | 51 kB | 26 | 2.6 kB | 159 | 159 | 2.5 kB | 377.00 us | - | +| | speed | 2.4 s | 65 kB | 49 kB | 35 | 3.6 kB | 109 | 126 | 1.7 kB | 343.00 us | +9% | +| | slots | 2.4 s | 66 kB | 53 kB | 29 | 2.9 kB | 106 | 159 | 1.7 kB | 391.00 us | -4% | +| g g > t t~ 3g | off | 6.2 s | 576 kB | 479 kB | 121 | 12.3 kB | 1890 | 1890 | 29.5 kB | 9.63 ms | - | +| | speed | 4.9 s | 463 kB | 408 kB | 213 | 21.6 kB | 1159 | 1551 | 18.1 kB | 8.97 ms | +7% | +| | slots | 5.0 s | 493 kB | 466 kB | 141 | 14.3 kB | 946 | 1890 | 14.8 kB | 9.93 ms | -3% | + +Two things worth reading off the amplitude columns. + +**The AMP array is recycled too, and it is where `slots` wins.** It used to be +declared at the full diagram count in every mode -- 113 kB at seven gluons, +with `speed` leaving 432 entries written by nobody -- which is what prompted +"Recycling the AMP array" below. Now `slots` runs 7245 amplitude calls through +946 entries at seven gluons, 14.8 kB rather than 113.2, while `speed` only +reaches 5869 because of the order it emits them in. + +**`slots` mode computes every amplitude** -- `amps computed` equals `amps +decl` on all eight of its rows. With no current sums nothing is skipped, so it +does the same amplitude work as the baseline *plus* the folds, and buys only a +shorter JAMP block and a shorter W array. That is why it is slower than off +almost everywhere rather than a wash: strictly more arithmetic for less +memory. `speed` is the opposite, skipping amplitudes outright (450 of 510, +6813 of 7245, 1551 of 1890), which is where its 5-8% comes from, and paying in +slots -- 121 to 213 at `g g > t t~ 3g`. + +**What it is actually good for.** Generation time and code size, more than +speed. `g g > 5 g` generates in 16.6 s rather than 35.2 s, a 53% cut and +reproducible: 385 seed diagrams unrolled is cheaper than 2485 generated. Its +`matrix.o` goes 5.8 MB to 3.4 MB and its `matrix.f` 3.3 MB to 2.5 MB. Runtime +is 4-9% above six particles and nothing at all below, and the `t t~` series +gains more than the pure gluon one at equal particle count, +9% at +`t t~ 2g` against +4% at `4g`. Peak RSS is flat except at seven gluons, +because the wavefunction store is a stack frame and the code image dominates. + +**How far it generalises.** Both series above are gluon-rich, so a third +process was measured as a check: `u u~ > z g g g g`, seven legs like +`g g > t t~ 3g` but with a Z and a quark line, so most of its diagrams have no +four gluon vertex at all. (`q` is not a defined multiparticle, hence `u u~`.) + +| process | mode | generate | matrix.f | matrix.o | W slots | AMP entries | amps computed | per call | speed | total/call | +|---|---|---|---|---|---|---|---|---|---|---| +| `g g > t t~ 3g` | off | 6.2 s | 576 kB | 479 kB | 121 | 1890 | 1890 | 9.63 ms | - | 41.8 kB | +| | speed | 4.9 s | 463 kB | 408 kB | 213 | 1159 | 1551 | 8.97 ms | +7% | 39.7 kB (-5%) | +| | slots | 5.0 s | 493 kB | 466 kB | 141 | 946 | 1890 | 9.93 ms | -3% | **29.1 kB (-30%)** | +| `u u~ > z 4g` | off | 2.9 s | 135 kB | 115 kB | 76 | 516 | 516 | 1.58 ms | - | 15.8 kB | +| | speed | 2.3 s | 127 kB | 110 kB | 85 | 391 | 450 | 1.51 ms | +4% | 14.7 kB (-7%) | +| | slots | 2.3 s | 132 kB | 117 kB | 84 | 384 | 516 | 1.65 ms | -5% | 14.5 kB (-8%) | + +**The payoff tracks the quartic fraction, exactly.** In `slots` mode the AMP +entry count is `total amplitudes - merge sources`, to within one, on every +process measured: + +| | amplitudes | merge sources | AMP in `slots` | saving | +|---|---|---|---|---| +| `g g > 5 g` | 7245 | 6300 (87%) | 946 | -87% | +| `g g > g g g g` | 510 | 405 (79%) | 106 | -79% | +| `g g > t t~ 3g` | 1890 | 945 (50%) | 946 | -50% | +| `u u~ > z 4g` | 516 | 132 (**26%**) | 384 | -26% | + +So `u u~ > z 4g` is the weakest case measured, and predictably: there is +simply little to merge. Runtime follows at +4% rather than +7 to +9%, and the +working set at -8% rather than -30 or -75%. + +It is also the one process where **`slots` does not reduce the wavefunctions +either** -- 84 against off's 76, worse -- which breaks the pattern from the +gluon-rich processes, where reversing the order always recovered them. On this +topology `slots` is the worst of the three: 5% slower than off and larger in +W, for an AMP count it barely wins over `speed`, 384 against 391. `speed` +still behaves, +4% with generation down 21% and a smaller source and object, +which is also what the `auto` gate picks at seven legs. + + +## The multiplicity gate + +The sweep says the full merging turns over at six external legs, and the same +sweep says it costs slots below that. So `auto` gates on it: +`Amplitude.generate_diagrams` only takes the seed rule when the process has +`madgraph.merge_quartic_min_legs` legs or more, six by measurement. `speed` +and `slots` asked for by name are unconditional -- that is how you get the +merging on a small process anyway. + +What is left below the threshold is not nothing, and this was worth measuring +rather than assuming. The *amplitude* merges do not need the seed rule: they +are found from the colour algebra by `unroll_quartic_vertices`, so they still +apply, and they shrink the JAMP block without touching the wavefunctions: + +| `g g > g g g` | slots | JAMP temporaries | per call | +|---|---|---|---| +| off | 12 | 72 | 87.75 us | +| `auto` (merges only) | 12 | 42 | **84.25 us**, +4.0% | +| `speed` (full) | 19 | 42 | 87.25 us, -1% | + +So below the threshold `auto` is *better* than both -- it keeps the JAMP fold, +which is free, and drops the reordering, which is what costs the seven extra +slots. At four and five legs elsewhere it is neutral rather than positive +(`g g > g g` 5.47 -> 5.48 us, `g g > t t~ g` 30.33 -> 30.42 us, both inside +the 0.6% floor) and never negative, at an unchanged slot count. + +Above the threshold nothing changes: `auto` at six legs generates a +byte-identical `matrix.f` to `speed`. + +## Recycling the AMP array + +`reuse_outdated_wavefunctions` recycles the wavefunctions; the amplitudes were +not recycled at all. `AMP` was declared `COMPLEX*16 AMP(NGRAPHS)` at the full +diagram count in every mode, so seven gluons allocated 113 kB of it and +`speed` left 432 entries written by nobody. + +`HelasMatrixElement.get_amplitude_slots` now does for AMP what +`reuse_outdated_wavefunctions` does for W. The enabling change is *where the +merges are written*: they used to be emitted in one block at the very end, +which kept every source alive to the end, and each is now written as soon as +both of its amplitudes exist. Once `AMP(t) = AMP(t) + AMP(s)` has run, `s` is +free. + +| process | mode | AMP entries | AMP | W slots | W | total per call | +|---|---|---|---|---|---|---| +| `g g > g g g` | off | 45 | 0.7 kB | 12 | 1.2 kB | 1.9 kB | +| | auto | 16 | 0.2 kB | 12 | 1.2 kB | 1.5 kB (-24%) | +| | speed | 24 | 0.4 kB | 19 | 1.9 kB | 2.3 kB (+20%) | +| | slots | 16 | 0.2 kB | 12 | 1.2 kB | 1.5 kB (-24%) | +| `g g > g g g g` | off | 510 | 8.0 kB | 51 | 5.2 kB | 13.1 kB | +| | speed | 316 | 4.9 kB | 78 | 7.9 kB | 12.9 kB (-2%) | +| | slots | **106** | 1.7 kB | 54 | 5.5 kB | 7.1 kB (**-46%**) | +| `g g > t t~ g g` | off | 159 | 2.5 kB | 26 | 2.6 kB | 5.1 kB | +| | speed | 109 | 1.7 kB | 35 | 3.6 kB | 5.3 kB (+3%) | +| | slots | 106 | 1.7 kB | 29 | 2.9 kB | 4.6 kB (-10%) | +| `g g > 5 g` | off | 7245 | 113.2 kB | 268 | 27.2 kB | 140.4 kB | +| | speed | 5869 | 91.7 kB | 259 | 26.3 kB | 118.0 kB (-16%) | +| | slots | **946** | 14.8 kB | 199 | 20.2 kB | **35.0 kB (-75%)** | + +**`slots` reaches the floor and `speed` does not**, and the reason is the +diagram order rather than anything about the allocator. Only (2n-5)!! of the +amplitudes are read by the JAMPs -- 105 at six gluons, 945 at seven -- and +everything else is a merge source which could in principle share a handful of +entries. `speed` emits every seed before its unrollings, so a source is born +early and its target arrives late and the entry cannot be reclaimed in +between: 5869 rather than 945. Reversing that order puts each source next to +its target, so `slots` lands on 946 and 106, one above the floor. + +That changes the `slots` case rather a lot. It used to buy 23% of the +wavefunction store and cost 6% more arithmetic; it now buys **75% of the whole +per-call working set** at seven gluons, which is the number that matters on a +gpu, where this is per thread. + +**It buys no time.** Measured at `g g > g g g g`, minimum of five: `speed` +2247 -> 2260 us and `slots` 2347 -> 2373 us across the change, both inside the +0.6% floor and if anything marginally the wrong way -- the arrays were already +cache resident at this size, and the merges are now interleaved rather than +batched. This is a memory optimisation, not a speed one. + +`NGRAPHS` only ever dimensioned `AMP` inside `matrix.f`, so it simply becomes +the entry count; `ngraphs.inc` keeps the diagram count. Everything reading AMP +afterwards goes through the same map: the JAMPs through +`ProcessExporterFortran.map_color_amplitudes`, and AMP2 through +`get_amplitude_slot_map`. AMP2 was the one to check, since multichannel reads +individual amplitudes -- it reads only the merge *targets*, which are the +entries that stay put, so it is unaffected. Verified on a madevent output at +six gluons: 316 entries written, 316 read, none read that is never written. + +`|M|^2` is unchanged on every row. + +## The diagram order — measured, and worth a lot + +`reuse_outdated_wavefunctions` is a linear scan allocator over lifetimes taken +in **emission order**, so `NWAVEFUNCS` depends on the order the diagrams are +written in. (The wavefunction *set* does not — that is content addressed, +`wavefunctions[wavefunctions.index(new_wf)]` — but the number of slots very +much does.) The shipped expansion emits every seed first and then the +unrollings, which is the worst case for it: a quartic current is made early +and its sum is not consumed until much later. + +Four orders were built and run. A sum can only be formed when its quartic +current is emitted before the target amplitude, which is what makes this a +trade rather than a free win: + +| `NWAVEFUNCS` / sums | off | seeds first (shipped) | by quartic count | seed then its own unrollings | last discovery | +|---|---|---|---|---|---| +| `g g > g g g` | 12 | 19 / 7 | 19 / 7 | 11 / 3 | **15 / 7** | +| `g g > g g g g` | 51 | 66 / 30 | 64 / 30 | 33 / **0** | **27 / 30** | +| `g g > 5 g` | 268 | 290 / 60 | 314 / 60 | 245 / 60 | **219 / 60** | + +Emitting a seed followed by its own unrollings gives the locality but loses +the sums: the fully cubic target is usually claimed by an earlier seed, so the +quartic current arrives after it. **Placing each diagram at its *last* +discovery instead of its first fixes that** — the target then sits after every +seed which can reach it — and takes both: all the sums, and a slot count +*below* the flag off baseline (27 against 51 at six gluons, 219 against 268 at +seven). + +**Shipped, once the reason it broke madmatrix was found.** The blocker was not +the order at all: + +1. The reconstruction could build the same cubic current with its two mothers + in either order. `sorted_mothers` leaves them alone, because for two + identical gluons its key ties and the sort is stable, so they stay two + objects. +2. **`VVV1P0_1` is antisymmetric under exchanging its two inputs.** Measured: + `VVV1P0_1(a,b) + VVV1P0_1(b,a) = 0` exactly. So the two are *negatives* of + each other and write out calls differing by a sign. +3. **`HelasWavefunction.__eq__` compares mothers by sorted number** and calls + them equal -- "the number for this wavefunction, the pdg code, and the + interaction id are irrelevant". +4. `export_cpp` renumbers wavefunctions by that equality, so the two landed on + one number and one slot inside a single matrix element and one of them + silently carried the wrong sign. The Fortran writer never hits it because + it does not renumber by equality. + +Taking the legs of both unrolled vertices in a canonical order removes every +such pair -- 20 of the 35 extra wavefunctions at six gluons, 30 of 57 at seven +-- so there is nothing left to collide, and the order goes in. + +| | wavefunctions | order-flipped twins | +|---|---|---| +| `g g > g g g` | 33 -> 39 | 0 | +| `g g > g g g g` | 111 -> 126 (was 146) | 20 -> 0 | +| `g g > 5 g` | 898 -> 925 (was 955) | 30 -> 0 | + +Placing each diagram at its last discovery then puts it after every seed which +can reach it, hence after every quartic current summable into it, so all the +sums survive. `NWAVEFUNCS` at seven gluons ends up **below** the unoptimised +build: 259 against 268. + +## The slot ordering, searched + +Once the twins are gone the slot count is the same in both backends (78 and +259 at six and seven gluons), so the order is one shared problem rather than a +per-backend one. Six orders were built and measured. The wavefunction count is +the same in all of them -- reordering the diagrams never changes *which* +currents exist, only how long each stays alive: + +| order | `g g > g g g` | `g g > g g g g` | `g g > 5 g` | sums kept | +|---|---|---|---|---| +| last discovery (shipped) | 19 | 78 | 259 | yes | +| first discovery | 19 | 81 | 314 | yes | +| by quartic count | 19 | 81 | 314 | yes | +| by quartic count, then last | 19 | 81 | 259 | yes | +| last discovery, then quartic count | 19 | 78 | 259 | yes | +| **reversed last discovery** | **12** | **54** | **199** | **no -- all lost** | + +The last row is the interesting one: it is far the best on slots and useless, +because reversing puts each target ahead of the quartic currents which feed +it, so no sum can be built. That is the trade in one line -- the constraint +that makes the sums possible is what costs the slots. + +A proper register-pressure greedy was also written: order the diagrams under +the precedence "everything which unrolls to a diagram comes before it", and at +each step take the one leaving fewest currents alive. It buys 19 -> 18 and +78 -> 76 and **nothing at all** at seven gluons, for an O(n^2) pass which +takes generation from 0.88 s to 3.24 s at seven gluons and would cost minutes +at eight. Not worth it; the shipped order is within a couple of slots of what +the search finds. + +## The backend-chosen order (`auto`) + +`speed` suits a cpu and `slots` suits a gpu, but one generation can feed both, +so `auto` defers the choice to `output`. It works only because the two modes +differ in nothing but the diagram order, and `slots` is `speed` reversed -- +verified byte for byte: generating in the `speed` order and reversing at +output time reproduces a native `slots` generation exactly, in both backends, +and a session going standalone -> standalone_mg7 -> standalone reproduces its +first output for the third. + +**The choice comes from the matrix element exporter, not the output format.** +`output madevent --me_exporter=` hands one +`self._curr_matrix_elements` to both exporters in a single `export_processes` +call, so the fortran driver and the gpu matrix elements *cannot* carry +different orders. Keying on the format would give that output the cpu order, +which is the one case the whole thing exists for; keying on the me exporter +gives it `slots` (measured: `NWAVEFUNCS=54` rather than 78 at six gluons). + +Two things had to be worked around, both pre-existing and neither specific to +this option: + +* **An export mutates the diagrams it is given** -- 345 of 757 vertex leg + records change on `g g > g g g g`. Reversing mutated diagrams gives an + equivalent but differently numbered result, so `apply_quartic_diagram_order` + reverses a copy taken before the first export. `copy.deepcopy` of the + *amplitudes* is not an option: it drags the model along and trips + `assert type(col_obj) != array.array` in `color_algebra.create_copy`. The + diagrams alone copy cleanly and cost 0.011 s at seven gluons. +* **An export also drops the marks saying the diagrams came from a seed** -- + after it, `seed_forbidden_cubic_ids` is empty and `quartic_unroll_tags` has + 0 entries instead of 405. So whether an amplitude may be reordered has to be + read before the first export, not at the output which wants to reorder. + +Two dead ends, both measured: reversing the `HelasMatrixElement` diagrams +instead of the amplitude's trips the lifetime assert in +`reuse_outdated_wavefunctions`, since a wavefunction number has to be first +seen in emission order; and restoring only `from_group` on the base diagrams +is not enough to undo an export. + +## Why this is not the default + +Defaulting `merge_quartic_vertices` to `auto` was tried and **reverted**. It is +safe on the paths it was built for -- UFO Fortran, madmatrix, the python +exporter -- and each of those is validated on `|M|^2`. Turning it on for +everything found three consumers which read the diagram or amplitude +*structure* rather than the result: + +1. **FKS born/real linking.** `link_rb_configs` finds the vertex splitting + `ij` into `i` and `j` and takes it out. The unrolling re-roots the real + diagrams and can put that pair in the closing vertex, where there is + nothing to take out, so the remainder is malformed and no born + configuration matches it. Same 3 diagrams selected, same tag set, different + decomposition: + + ``` + off ((1,2>1),(4,5>4),(1,3,4)) the 4-5 vertex is internal + auto ((1,2>1),(1,3>1),(1,4,5)) the 4-5 pair closes the diagram + ``` + + `p p > j j [QCD]` raised `FKSProcessError`. **Fixed** by generating an NLO + process with the merging off, in `FKSMultiProcess.__init__`, so the option + is now safe for an NLO user rather than only for the default. `g g > g g + [QCD]` generates byte-identically with the option set and unset. + +2. **The legacy `FortranHelasCallWriter`.** Only `FortranUFOHelasCallWriter` + emits the amplitude folds which put the merged contributions back. The + MG4-style writer computes `AMP(1..3)` from `GGGGXX` and then leaves them out + of the JAMPs -- a **silently wrong** `|M|^2`, not a crash. Not fixed: like + `export_cpp` and `export_python` it needs `merge_quartic_amplitudes=False`, + but `get_JAMP_lines` is on the exporter and does not know which writer it + is paired with. + +3. **Anything pinning the diagram order.** Cosmetic but wide: `colorize` and + `DiagramTag` tests select diagrams by position, and the sextet colour basis + goes 13 -> 15 because folding amplitudes decomposes the same `|M|^2` over + more colour structures (`|M|^2` bit-identical, checked with + `MatrixElementEvaluator`). + +### The scan for a fourth + +Looked for one, did not find one in the shipped tree, and the search bounds +the remaining risk. + +**Every merged-JAMP consumer, against every writer.** `get_color_amplitudes` +has four call sites: `export_cpp` (x2), `export_python` and `madmatrix` all +pass `merge_quartic_amplitudes=False`; only `export_v4`'s three +`get_JAMP_lines*` take the merged default, and every Fortran exporter pairs +with `FortranUFOHelasCallWriter`, which emits the folds. The base +`get_amplitude_merge_lines` returns `[]` and `FortranHelasCallWriter` does not +override `get_matrix_element_calls`, so it is the one writer that silently +drops them -- and it is selected exactly when `self._model_v4_path` is set, +i.e. under `import model_v4`. No other combination in the tree reaches merged +JAMPs without folds. + +The property which keeps `merge_quartic_amplitudes=False` safe is that +`get_color_amplitudes` drops the current-sum folded amplitudes +*unconditionally* and only the amplitude merges conditionally -- so a writer +which emits the sums but not the folds still gets consistent JAMPs. + +**A mechanical audit of the generated code.** For each `AMP(n)` in a generated +matrix element, whether it is written and whether it is read. Read-never- +written is garbage; written-never-read is a contribution dropped on the floor, +which is the legacy writer's signature. Run with the flag off as a control on +every output -- standalone, matchbox, madevent grouped and not, a decay chain, +helicity-recycled files, `u u~ > g g g`, `g g > t t~ g g`, `u u~ > u u~ g g`, +four to six gluons -- and clean everywhere. The split order path is included: +`ProcessExporterFortranSA` and `ProcessExporterFortranME` both always go +through `get_JAMP_lines_split_order`, whose `amp_orders` lists folded +amplitude numbers, but they never reach the code because the colour amplitudes +no longer mention them. + +**The two places which weight or group results.** Subprocess grouping for +`p p > j j` gives the same five directories and byte-identical `configs.inc` +and `coloramps.inc`; only the `g g > g g` matrix element differs. +`find_symmetry`, which feeds the multiplicative `symfact.dat`, keeps the same +equivalence classes and multiplicities -- `[3, 3, 3, 6]` at five gluons and +`[3, 6, 12, 12, 12, 12, 12, 12, 24]` at six, in both settings, with the same +number of channels and the same total weight. Only the representative diagram +indices renumber, which is the renumbering itself. + +So the residual risk of this class is two named things: `import model_v4`, and +a third-party plugin supplying its own `helas_exporter` paired with +`export_v4`'s merged JAMPs. + +The pattern is that the optimisation changes the *representation* -- diagram +order, rooting, which amplitudes survive into the JAMPs -- and every consumer +which reads representation rather than result has to be checked. Three turned +up in one pass, so defaulting it on wants an audit of those consumers, not +another round of patching outward. + +Also fixed on the way, and independent of all this: `link_rb_configs` built +`real_tags` deduplicated but left `good_diags` as it was, then walked the two +in lockstep -- `real_tags.remove(btag)` beside `good_diags.pop(ir)` -- so they +only stayed aligned while the dedup dropped nothing. A no-op on every process +in the test suite, but it made the result order dependent for no reason. + +## Where to go next + +**Give madmatrix the amplitude sums too.** It is the only backend without +them, because there is no `AMP` array to fold into — each amplitude goes +straight into the JAMPs. That is why it gains 1.9% where Fortran gains 3.9%, +and why five gluons still loses. The sources for one target could be +accumulated into a second `amp_sv` slot before the JAMP lines are written, +which the seed ordering makes possible (the quartic diagrams come first), at +the cost of one accumulator per open target. + +Then there are the sums which do not sit at an amplitude, and they need a +node to have exactly one rooting *per merge*, which a diagram list cannot +give. Three ways on, in increasing size: + +1. **The second substitution.** Cheapest of the three, ceiling 432 amplitude + calls at seven gluons, but measured to be at most 282 of them and possibly + much less. Two things stand in the way and only the first is bookkeeping: + + *Identification.* All 432 targets which have a two-substitution source have + exactly one of the two singles present; the other is the same diagram + rooted differently, so it is a different amplitude object with different + mothers and `match_quartic_mothers` cannot see it. Finding it by its + diagram and colour chain — the frame `compute_quartic_amplitude_merges` + already works in — rather than by its mothers would take it. + + *The coefficients have to multiply.* Substituting two mothers also produces + the amplitude with both substituted, weighing the product of the two + coefficients, so the merge map has to agree. It does not always: taking the + double's coefficient over the known single's, the missing single would have + to weigh -1 for 360 of the 432 targets and +1 for 72, and for **150 of them + no merge source into that target weighs that at all**. The sign from + `diagram_colour_signature` does not factorise over two contractions in + general. `subset_is_merged` already refuses those, which is what keeps the + present code sound; extending the identification does not remove the check. +2. **Partial CSE.** Keep the DAG, add `TMP = W1 + W4` alongside W1 and W4, + and split the consumers. Bounded gain: at six gluons only 2 of the 6 + quartic consumers correspond, so it saves 2 subtrees per node out of 50 + nodes. +3. **Drop the diagram list for the currents.** Build the wavefunctions by a + Berends-Giele recursion over subsets — at each node, cubic pair plus + quartic, which generates exactly the matchings and never double counts — + and keep the 220 diagrams only for what they are actually needed for + (multichannel, `matrix.ps`). + +Anything that fragments the diagram list to get a per-rooting copy runs into +pitfall 1, and anything that expands each seed independently double counts — +the fully cubic diagrams get reached once per matching of their adjacency +graph (225 instead of 105 at six gluons). + +## Pitfalls — all of these cost real time in the previous session + +1. **Do not fragment the diagram list.** Generating one diagram per colour + structure (220 -> 510) works numerically but changes a user-visible number + and the MadEvent multichannel. It was reverted for that reason. +2. **`copy.deepcopy` on a `Model` silently breaks colour.** `ColorObject` + derives from `array.array` and deepcopy degrades it to a plain `array`, so + the structures stop being recognisable. Use `create_copy`. +3. **Two colour-chain conventions exist.** `helas_objects` colorizes the + *reconstructed* amplitude (`get_base_amplitude`), whose vertex order can + differ from the generated diagrams. They agree up to six gluons and diverge + at seven. Anything mapping chains to amplitudes must use the reconstructed + one. +4. **Pinning a single colour structure needs colour and Lorentz in the same + frame.** The sum over the three structures is permutation invariant but the + individual terms are not, so the structure must be chosen against + `sorted_mothers` (what ALOHA receives), not the vertex leg order. +5. **The `get_color_amplitudes` filter and the writer emitting the folds must + land together.** Filtering without emitting silently drops the quartic + contributions from the JAMPs — this produced a wrong `|M|^2` that survived + three rounds of debugging because the helas calls looked byte-identical to + baseline. When a number is wrong, diff the JAMP block first. +6. **A current sum is only legitimate if the consumers correspond 1:1.** + Summing into a shared current hands the contribution to *every* consumer. + Measured counter-example on the fragmented structure: quartic current 11 + had consumers {7,9,36,38,60,62} mapping to {1,31,55}, while its partner had + {1,5,31,34,55,58} — three consumers would have gained a term they must not + have. +7. **The amplitude sums buy nothing at runtime.** The JAMP CSE already finds + those pairs (`TMP_JAMP(2) = AMP(1) + AMP(4)`). They do shrink the JAMP block + (1091 -> 697 lines at six gluons), which helps the optimiser and compile + time. Expect the speedup to come from the currents, not from these. + *Measured since:* wrong at seven gluons, where the JAMP block goes 23672 -> + 7864 lines and that alone is +3.1%. Right at six, where it is a wash. +8. **A wavefunction number is not a slot.** `reuse_outdated_wavefunctions` + hands the same `me_id` to wavefunctions whose lifetimes do not overlap, and + works those lifetimes out from the diagrams alone. Anything emitting an + extra read of a wavefunction has to extend the lifetime there, or it reads + whatever else has since been written into that slot -- which showed up as + `W(11)%W(:) = W(11)%W(:) + W(11)%W(:)`. + +## Measuring + +Generation: `set merge_quartic_vertices speed` then `generate g g > g g g g`. +Compare `matrix.f` against a run without the variable. Tests: +`./tests/test_manager.py -p U test_diagram_generation test_color_amp +test_helas_objects test_base_objects` (199, must stay green with the flag off). diff --git a/madgraph/__init__.py b/madgraph/__init__.py index 0d17042f0..000a1794c 100755 --- a/madgraph/__init__.py +++ b/madgraph/__init__.py @@ -60,4 +60,42 @@ class aMCatNLOError(MadGraph5Error): ordering = True else: ordering = False + +# Sum the quartic gluon contributions into the cubic amplitude carrying the +# same colour factor, see HelasMatrixElement.get_quartic_amplitude_merges. +# Set through the interface, "set merge_quartic_vertices ", and read +# here because it is wanted while the diagrams are generated, long before any +# exporter exists. False, or one of: +# 'speed' -- the current sums, and the diagram order which allows them. Wins +# on cpu, where the amplitude calls dominate. +# 'slots' -- no current sums, and the order which keeps fewest currents +# alive. Trades 6% more amplitude calls for 23% fewer +# wavefunction slots at seven gluons, which is the trade a gpu +# wants when occupancy is the limit. +# 'auto' -- decide per process. Below merge_quartic_min_legs external legs +# only the amplitude merges are taken, which shrink the JAMP +# block for free; the seed rule -- and with it the reordering, +# the current sums and the slots they cost -- is left off, since +# below that size it does not pay for itself. At or above it, +# generate as for 'speed' and let each output pick 'slots' when +# the matrix elements go to a gpu backend and 'speed' otherwise. +# The interface resolves the second half before anything reads it +# again, so only the generation ever sees 'auto'. +# 'speed' and 'slots' are unconditional -- they are the way to ask for the +# merging on a small process anyway. +# Off by default: several consumers read the diagram or amplitude structure +# rather than the result, see "Why this is not the default" in +# docs/gluon-quartic-plan.md. +merge_quartic_vertices = False + +# External legs an 'auto' process needs before the seed rule is worth it. +# Measured on g g > N g and g g > t t~ N g, both of which turn over at six: +# five legs and below the full merging is inside the noise or slightly +# negative (-1% at g g > g g g, which also pays 12 -> 19 wavefunction slots), +# six and above it is a steady 5-8% with the source and the object file +# shrinking too. The amplitude merges alone, which is what is left below the +# threshold, are +4% at g g > g g g and neutral at four and five legs +# elsewhere, at the same slot count -- never a loss. See "Full sweep" in +# docs/gluon-quartic-plan.md. +merge_quartic_min_legs = 6 diff --git a/madgraph/core/base_objects.py b/madgraph/core/base_objects.py index a43e58fad..18eff6a7b 100755 --- a/madgraph/core/base_objects.py +++ b/madgraph/core/base_objects.py @@ -2980,7 +2980,7 @@ class Vertex(PhysicsObject): """ sorted_keys = ['id', 'legs'] - + # This sets what are the ID's of the vertices that must be ignored for the # purpose of the multi-channeling. 0 and -1 are ID's of various technical # vertices which have no relevance from the perspective of the diagram @@ -3014,7 +3014,6 @@ def default_setup(self): # that it can be easily identified when constructing the DiagramChainLinks. self['id'] = 0 self['legs'] = LegList() - def filter(self, name, value): """Filter for valid vertex property values.""" diff --git a/madgraph/core/color_amp.py b/madgraph/core/color_amp.py index 55d5731fe..20850bde9 100755 --- a/madgraph/core/color_amp.py +++ b/madgraph/core/color_amp.py @@ -197,7 +197,7 @@ def add_vertex(self, vertex, diagram, model, new_res_dict = {} for i, col_str in \ enumerate(inter_color): - + # Ignore color string if it doesn't correspond to any coupling if i not in inter_indices: continue diff --git a/madgraph/core/diagram_generation.py b/madgraph/core/diagram_generation.py index 1f239abf5..f5ce93b80 100755 --- a/madgraph/core/diagram_generation.py +++ b/madgraph/core/diagram_generation.py @@ -24,11 +24,13 @@ import array import copy +import fractions import itertools import logging import madgraph import madgraph.core.base_objects as base_objects +import madgraph.core.color_algebra as color_algebra import madgraph.various.misc as misc import madgraph.fks.fks_tag as fks_tag from madgraph import InvalidCmd, MadGraph5Error @@ -427,6 +429,435 @@ def __str__(self): __repr__ = __str__ +#=============================================================================== +# Unrolling of quartic vertices into pairs of cubic vertices +#=============================================================================== + +class UnrollDiagramTag(DiagramTag): + """DiagramTag which keeps every external leg distinct. + + The default DiagramTag deliberately identifies identical final state + particles, which is what is wanted when comparing different processes. To + recognise a diagram inside its own amplitude we need the leg numbers, so + that e.g. the t- and u-channel of g g > g g do not collide.""" + + @staticmethod + def link_from_leg(leg, model): + return [((leg.get('id'), leg.get('number')), leg.get('number'))] + + +def get_unrollable_quartic_vertices(model): + """Find the quartic interactions which factorise into two cubic ones. + + A quartic vertex qualifies when its four legs carry the same + self-conjugate particle and each of its colour structures is a product of + two structure constants sharing a single summed index -- the four gluon + vertex being the canonical example. Such a colour structure splits the + four colour indices into two pairs, and that splitting is precisely the + one produced by two cubic vertices joined by an internal line. + + Returns {quartic_id: (cubic_id, pairings)} where pairings[icolor] is the + pair of 2-tuples of colour index positions separated by the summed index. + """ + + try: + return model._unrollable_quartic_vertices + except AttributeError: + pass + + # Cubic candidates: three identical legs with a single f(0,1,2) structure + cubic_by_pdg = {} + for inter in model.get('interactions'): + parts = inter.get('particles') + if len(parts) != 3: + continue + pdgs = misc.make_unique([p.get_pdg_code() for p in parts]) + if len(pdgs) != 1: + continue + color = inter.get('color') + if len(color) != 1 or len(color[0]) != 1: + continue + if not isinstance(color[0][0], color_algebra.f) or \ + sorted(color[0][0]) != [0, 1, 2]: + continue + cubic_by_pdg[pdgs[0]] = inter.get('id') + + res = {} + for inter in model.get('interactions'): + parts = inter.get('particles') + if len(parts) != 4: + continue + pdgs = misc.make_unique([p.get_pdg_code() for p in parts]) + if len(pdgs) != 1 or pdgs[0] not in cubic_by_pdg: + continue + pairings = [] + for col_str in inter.get('color'): + pairing = _f_pair_split(col_str) + if pairing is None: + break + pairings.append(pairing) + else: + if pairings: + res[inter.get('id')] = (cubic_by_pdg[pdgs[0]], pairings) + + try: + model._unrollable_quartic_vertices = res + except AttributeError: + pass + return res + + +def _f_pair_split(col_str): + """If col_str is f(..)*f(..) with a single shared summed index and the + four remaining indices being 0,1,2,3, return the two pairs of external + colour index positions it separates, otherwise None.""" + + if len(col_str) != 2 or not all(isinstance(obj, color_algebra.f) + for obj in col_str): + return None + shared = [i for i in col_str[0] if i in col_str[1]] + if len(shared) != 1 or shared[0] >= 0: + return None + pairs = tuple(tuple(sorted(i for i in obj if i != shared[0])) + for obj in col_str) + if sorted(pairs[0] + pairs[1]) != [0, 1, 2, 3]: + return None + return pairs + + +def glued_vertices(diagram): + """Vertices of the diagram with the trailing identity vertex glued into + the one before it. + + The identity vertex only states that its two legs are the two ends of one + and the same line, so it is dropped once the diagram is complete, by + handing that line to the vertex before it. Until that is done the same + diagram can be written in several ways, and comparing two of them is + meaningless. Returns None when there is nothing to glue.""" + + vertices = diagram.get('vertices') + if len(vertices) <= 1 or vertices[-1].get('id') != 0: + return None + + vertices = copy.copy(vertices) + lastvx = vertices.pop() + nexttolastvertex = copy.copy(vertices.pop()) + legs = copy.copy(nexttolastvertex.get('legs')) + ntlnumber = legs[-1].get('number') + lastleg = [leg for leg in lastvx.get('legs') + if leg.get('number') != ntlnumber][0] + # Reset onshell in case we have forbidden s-channels + if lastleg.get('onshell') == False: + lastleg.set('onshell', None) + # Replace the last leg of nexttolastvertex + legs[-1] = lastleg + nexttolastvertex.set('legs', legs) + vertices.append(nexttolastvertex) + return vertices + + +def get_unrollable_cubic_ids(model): + """Interaction ids of the cubic vertices a factorisable quartic vertex + unrolls into -- the three gluon vertex for the four gluon one.""" + + return frozenset(cubic_id for cubic_id, pairings in + get_unrollable_quartic_vertices(model).values()) + + +def colour_index_order(vertex, is_last, model): + """Return the vertex legs in the order in which color_amp.ColorBasis maps + them onto the colour indices of the interaction. + + This mirrors ColorBasis.add_vertex: the outgoing leg of an internal + vertex is flipped to its antiparticle and moved to the front, then the + legs are sorted following the particle order of the interaction. Note + that the sorting can move the outgoing leg away from the front again, so + its position is returned along with the legs. + + Returns (legs, outgoing_position), the position being None for the last + vertex of a diagram, or None if the vertex does not match its + interaction.""" + + inter = model.get_interaction(vertex.get('id')) + if inter is None: + return None + + legs = vertex.get('legs') + entries = [] + for index, leg in enumerate(legs): + part = model.get('particle_dict')[leg.get('id')] + outgoing = index == len(legs) - 1 and not is_last + entries.append((leg, part.get_anti_pdg_code() if outgoing + else part.get_pdg_code(), outgoing)) + + if not is_last: + entries.insert(0, entries.pop(-1)) + + ordered = [] + for pdg in [p.get_pdg_code() for p in inter.get('particles')]: + for index, entry in enumerate(entries): + if entry[1] == pdg: + ordered.append(entries.pop(index)) + break + else: + return None + + out_position = None + for index, entry in enumerate(ordered): + if entry[2]: + out_position = index + + return [entry[0] for entry in ordered], out_position + + +def split_quartic_vertex(vertex, is_last, pairing, cubic_id, model): + """Replace a quartic vertex by the two cubic vertices that the given + colour structure factorises into, joined by a new internal line.""" + + ordered, out_position = colour_index_order(vertex, is_last, model) + first, second = pairing + if not is_last and out_position in first: + # the pair the outgoing leg does not belong to is the one replaced by + # the new internal line + first, second = second, first + + combined = sorted([ordered[i] for i in first], + key=lambda leg: leg.get('number')) + new_leg = base_objects.Leg({ + 'id': combined[0].get('id'), + 'number': min(leg.get('number') for leg in combined), + 'state': len([l for l in combined if not l.get('state')]) != 1, + 'from_group': True}) + + first_vx = base_objects.Vertex({ + 'legs': base_objects.LegList(combined + [new_leg]), + 'id': cubic_id}) + if is_last: + rest = sorted([ordered[i] for i in second] + [new_leg], + key=lambda leg: leg.get('number')) + else: + # the outgoing leg has to stay last + rest = sorted([ordered[i] for i in second if i != out_position] + + [new_leg], key=lambda leg: leg.get('number')) + \ + [ordered[out_position]] + second_vx = base_objects.Vertex({'legs': base_objects.LegList(rest), + 'id': cubic_id}) + + return [first_vx, second_vx] + + +# Placeholder standing for the single summed index of a factorisable quartic +# colour structure, which is not a leg of the vertex. +_SUMMED = 'summed' + + +def reroot_diagram(diagram, root, model): + """Return the same diagram decomposed around another final vertex. + + A diagram is a tree, and which of its vertices is written last decides + which internal lines become currents: everything on the far side of the + final vertex is built up as a wavefunction, while the final vertex only + produces an amplitude. Re-rooting therefore changes which currents exist + without changing the diagram -- same vertices, same lines, same physics -- + which is what lets a quartic current find the cubic current it has to be + summed with. + + root is the index of the vertex to end on. Returns None if the diagram + does not decompose into a tree, which should not happen. + """ + + vertices = diagram.get('vertices') + last = len(vertices) - 1 + if root == last: + return diagram + + # Split every vertex into the external legs it holds and the internal + # lines it shares with another vertex. A line is recognised by its number + # being live, since a vertex reuses the smallest incoming number for the + # leg it produces. + live = {} + externals = [[] for _ in vertices] + lines = [] + for i, vertex in enumerate(vertices): + incoming = vertex.get('legs') if i == last else vertex.get('legs')[:-1] + for leg in incoming: + producer = live.pop(leg.get('number'), None) + if producer is None: + externals[i].append(leg) + else: + lines.append((producer, i, + vertices[producer].get('legs')[-1])) + if i != last: + live[vertex.get('legs')[-1].get('number')] = i + if live: + return None + + neighbours = {} + line_of = {} + for producer, consumer, leg in lines: + neighbours.setdefault(producer, []).append(consumer) + neighbours.setdefault(consumer, []).append(producer) + line_of[(producer, consumer)] = (leg, True) + line_of[(consumer, producer)] = (leg, False) + + # Walk out from the new root so that every vertex is emitted after the + # ones now feeding it. + order = [] + parent = {root: None} + def visit(node): + for other in neighbours.get(node, []): + if other == parent[node]: + continue + parent[other] = node + visit(other) + order.append(node) + visit(root) + if len(order) != len(vertices): + return None + + produced = {} + res = base_objects.VertexList() + for node in order: + incoming = list(externals[node]) + for other in neighbours.get(node, []): + if other != parent[node]: + incoming.append(produced[other]) + legs = base_objects.LegList(incoming) + if node != root: + leg, forwards = line_of[(node, parent[node])] + part = model.get('particle_dict')[leg.get('id')] + outgoing = base_objects.Leg({ + # the line keeps its particle if it already pointed this way, + # and is flipped when the re-rooting reverses it + 'id': leg.get('id') if forwards or part.get('self_antipart') + else -leg.get('id'), + 'number': min(l.get('number') for l in incoming), + 'state': len([l for l in incoming + if not l.get('state')]) != 1, + 'from_group': True}) + produced[node] = outgoing + legs.append(outgoing) + res.append(base_objects.Vertex({'legs': legs, + 'id': vertices[node].get('id')})) + + return base_objects.Diagram({'vertices': res}) + + +def diagram_colour_signature(diagram, model, color_chain, unrollable): + """Canonical signature of the colour string of one colour structure choice. + + Returns (key, coeff), where key identifies the product of colour objects + up to the antisymmetry of the structure constants and coeff collects the + resulting sign together with the rational prefactors. Two contributions + sharing a key have proportional colour factors and can be summed, the + relative weight being the ratio of their coeff. + + Every line is labelled by the set of external legs it separates rather + than by its leg number, so that the signature can be compared between + diagrams which number their internal lines differently. Returns None when + the colour structure is not of a supported form. + """ + + vertices = diagram.get('vertices') + last = len(vertices) - 1 + + # A leg number is reused by the vertex which produces it, so an external + # leg is one consumed while no internal line of that number is alive. + alive = {} + externals = [] + for i, vertex in enumerate(vertices): + legs = vertex.get('legs') + incoming = legs if i == last else legs[:-1] + for leg in incoming: + if alive.pop(leg.get('number'), None) is None: + externals.append(leg.get('number')) + if i != last: + alive[legs[-1].get('number')] = True + externals = frozenset(externals) + if not externals: + return None + reference = min(externals) + + def normalise(side): + """Root independent label of the line splitting side from its rest.""" + return side if reference not in side else externals - side + + factors = [] + coeff = fractions.Fraction(1, 1) + imaginary = False + live = {} + for i, vertex in enumerate(vertices): + inter = model.get_interaction(vertex.get('id')) + if inter is None: + return None + order = colour_index_order(vertex, i == last, model) + if order is None: + return None + ordered, out_position = order + sides = [None] * len(ordered) + below = frozenset() + # the outgoing leg is only known once all the incoming ones are read, + # since a vertex reuses the smallest incoming leg number for it + for pos in range(len(ordered)): + if pos == out_position: + continue + number = ordered[pos].get('number') + side = live.pop(number, None) + if side is None: + side = frozenset([number]) + sides[pos] = side + below = below | side + if out_position is not None: + # the outgoing leg carries everything that is not below it + sides[out_position] = externals - below + live[ordered[out_position].get('number')] = below + + labels = dict(enumerate(normalise(side) for side in sides)) + if vertex.get('id') in unrollable: + pairing = unrollable[vertex.get('id')][1][color_chain[i]] + labels[_SUMMED] = normalise(sides[pairing[0][0]] | + sides[pairing[0][1]]) + + if not inter.get('color'): + continue + col_str = inter.get('color')[color_chain[i]] + coeff *= col_str.coeff + imaginary ^= col_str.is_imaginary + for obj in col_str: + try: + indices = [labels[j if j >= 0 else _SUMMED] for j in obj] + except KeyError: + return None + if isinstance(obj, color_algebra.f): + # totally antisymmetric: sort and keep track of the parity + order = sorted(range(len(indices)), + key=lambda k: sorted(indices[k])) + coeff *= _permutation_sign(order) + indices = [indices[k] for k in order] + factors.append((obj.__class__.__name__, + tuple(tuple(sorted(index)) for index in indices))) + + return (tuple(sorted(factors)), imaginary), coeff + + +def _permutation_sign(order): + """Signature of a permutation given as a list of positions.""" + + sign = 1 + seen = [False] * len(order) + for start in range(len(order)): + if seen[start]: + continue + length = 0 + pos = start + while not seen[pos]: + seen[pos] = True + pos = order[pos] + length += 1 + if length % 2 == 0: + sign = -sign + return sign + #=============================================================================== # Amplitude #=============================================================================== @@ -436,6 +867,13 @@ class Amplitude(base_objects.PhysicsObject): generate the diagrams for the amplitude """ + # Interaction ids of the cubic vertices which the seed rule forbids to + # share a line, see generate_diagrams. Empty -- so the rule is inactive -- + # unless madgraph.merge_quartic_vertices is set. + seed_forbidden_cubic_ids = frozenset() + # Links recorded while the seed was expanded, see expand_seed_diagrams + quartic_unroll_tags = {} + def default_setup(self): """Default values for all properties""" @@ -580,8 +1018,7 @@ def generate_diagrams(self, returndiag=False, diagram_filter=False): "particles are missing in model: %s" % model.get('particles') assert model.get('interactions'), \ - "interactions are missing in model" - + "interactions are missing in model" res = base_objects.DiagramList() # First check that the number of fermions is even @@ -668,6 +1105,23 @@ def generate_diagrams(self, returndiag=False, diagram_filter=False): max_multi_to1 = max([len(key) for key in \ model.get('ref_dict_to1').keys()]) + # Seed rule: when the quartic vertices are to be put back afterwards + # by unrolling (see unroll_quartic_vertices), generate only the + # diagrams unrolling cannot produce. Unrolling a quartic vertex always + # yields two cubic vertices sharing the line which replaced it, so a + # diagram can be reconstructed exactly when two of its cubic vertices + # share a line -- and the seed is what is left over. + # Left off for a decay chain, whose identity vertex is kept rather + # than glued in, and for loop amplitudes, whose diagram set is not the + # one the unrolling reasons about. 'auto' also leaves it off below + # merge_quartic_min_legs, where it is measured to buy nothing. + self.seed_forbidden_cubic_ids = frozenset() + if madgraph.merge_quartic_vertices and not self.has_loop_process() \ + and not process.get('is_decay_chain') \ + and (madgraph.merge_quartic_vertices != 'auto' or + len(process.get('legs')) >= madgraph.merge_quartic_min_legs): + self.seed_forbidden_cubic_ids = get_unrollable_cubic_ids(model) + # Reduce the leg list and return the corresponding # list of vertices @@ -707,6 +1161,10 @@ def generate_diagrams(self, returndiag=False, diagram_filter=False): for vertex_list in reduced_leglist: res.append(self.create_diagram(base_objects.VertexList(vertex_list))) + # Put back the diagrams the seed rule left out + if self.seed_forbidden_cubic_ids: + res = self.expand_seed_diagrams(res) + # Record whether or not we failed generation before required # s-channel propagators are taken into account failed_crossing = not res @@ -822,25 +1280,11 @@ def generate_diagrams(self, returndiag=False, diagram_filter=False): # Replace final id=0 vertex if necessary if not process.get('is_decay_chain'): for diagram in res: - vertices = diagram.get('vertices') - if len(vertices) > 1 and vertices[-1].get('id') == 0: - # Need to "glue together" last and next-to-last - # vertex, by replacing the (incoming) last leg of the - # next-to-last vertex with the (outgoing) leg in the - # last vertex - vertices = copy.copy(vertices) - lastvx = vertices.pop() - nexttolastvertex = copy.copy(vertices.pop()) - legs = copy.copy(nexttolastvertex.get('legs')) - ntlnumber = legs[-1].get('number') - lastleg = [leg for leg in lastvx.get('legs') if leg.get('number') != ntlnumber][0] - # Reset onshell in case we have forbidden s-channels - if lastleg.get('onshell') == False: - lastleg.set('onshell', None) - # Replace the last leg of nexttolastvertex - legs[-1] = lastleg - nexttolastvertex.set('legs', legs) - vertices.append(nexttolastvertex) + # "glue together" last and next-to-last vertex, by replacing + # the (incoming) last leg of the next-to-last vertex with the + # (outgoing) leg in the last vertex + vertices = glued_vertices(diagram) + if vertices is not None: diagram.set('vertices', vertices) if res and not returndiag: @@ -944,9 +1388,245 @@ def remove_diag(diag, model=None): return res + def expand_seed_diagrams(self, seed): + """Put back the diagrams the seed rule left out of the generation. + + A seed diagram stands for itself and for every diagram obtained by + replacing any subset of its quartic vertices by the pair of cubic + vertices one of their colour structures factorises into. Several seeds + reach the same diagram, hence the dedup, and the result is the full + diagram set again -- same diagrams, same count, so nothing + user-visible moves. + + What the detour buys is the decomposition: an unrolled diagram is the + seed with one vertex taken apart, so it is rooted exactly like the + diagram it has to be summed with and shares every current except the + one being summed. + + Also records the link between a quartic contribution and the diagram + it unrolls to, see get_quartic_unroll_links. + """ + + model = self.get('process').get('model') + unrollable = get_unrollable_quartic_vertices(model) + ninitial = self.get_ninitial() + + def canonical_tag(diagram): + return str(UnrollDiagramTag(diagram, model, ninitial)) + + res = base_objects.DiagramList() + seen = set() + tag_of = [] + last_seen = {} + clock = [0] + + def touch(tag): + # A diagram is placed at its *last* discovery, so that it lands + # after every seed which can reach it -- and so after every + # quartic current which can be summed into it. Seeing it again + # only moves it later, it is never generated twice. + clock[0] += 1 + last_seen[tag] = clock[0] + + self.quartic_unroll_tags = {} + todo = [] + for diagram in seed: + # The identity vertex is glued in right away rather than at the + # end of generate_diagrams. Until it is, the same diagram has + # several spellings and comparing two of them is meaningless, and + # a quartic vertex sitting just before it is not yet the last one + # -- which is what decides how its colour structures are indexed. + vertices = glued_vertices(diagram) + if vertices is not None: + diagram = self.create_diagram(vertices) + tag = canonical_tag(diagram) + seen.add(tag) + touch(tag) + res.append(diagram) + tag_of.append(tag) + todo.append((diagram, tag)) + + # Unrolling is confluent, so taking the diagrams it produces through + # the same treatment adds nothing to the set -- but it does give the + # link for the quartic vertices they have left. + cursor = 0 + while cursor < len(todo): + diagram, own_tag = todo[cursor] + cursor += 1 + vertices = diagram.get('vertices') + positions = [i for i, vertex in enumerate(vertices) + if vertex.get('id') in unrollable] + if not positions: + continue + # None leaves the vertex alone, the other entries are the colour + # structures carrying a coupling, as kept by ColorBasis.colorize + allowed = [[None] + sorted(misc.make_unique( + [key[0] for key in model.get_interaction( + vertices[p].get('id')).get('couplings')])) + for p in positions] + + for keys in itertools.product(*allowed): + choice = dict((position, key) for position, key + in zip(positions, keys) if key is not None) + if not choice: + continue + unrolled = self.unrolled_diagram(diagram, choice, unrollable) + tag = canonical_tag(unrolled) + touch(tag) + if tag not in seen: + seen.add(tag) + res.append(unrolled) + tag_of.append(tag) + todo.append((unrolled, tag)) + if len(choice) == len(positions): + chain = tuple(choice.get(i, 0) + for i in range(len(vertices))) + self.quartic_unroll_tags[(own_tag, chain)] = tag + + order = sorted(range(len(res)), key=lambda i: (last_seen[tag_of[i]], i)) + if madgraph.merge_quartic_vertices == 'slots': + # Reversed, each diagram lands at its *first* discovery instead. + # That keeps far fewer currents alive -- 199 slots against 259 at + # seven gluons -- at the price of the current sums, since a target + # then comes before the quartic currents which feed it. The trade + # a gpu wants, where the wavefunction store is per thread. + order.reverse() + return base_objects.DiagramList([res[i] for i in order]) + + def get_quartic_unroll_links(self, diaglist=None): + """Return {(diagram index, colour chain): target index} recorded while + the seed was expanded, resolved against the diagram list as it stands. + + This is the same map as the one unroll_quartic_vertices reconstructs + from the colour algebra, but obtained for free: the two diagrams are + one and the same seed unrolled differently. Empty when the diagrams + were not generated from a seed.""" + + if not self.quartic_unroll_tags: + return {} + + model = self.get('process').get('model') + if diaglist is None: + diaglist = self.get('diagrams') + ninitial = self.get_ninitial() + index = dict((str(UnrollDiagramTag(diagram, model, ninitial)), i) + for i, diagram in enumerate(diaglist)) + + res = {} + for (seed_tag, chain), tag in self.quartic_unroll_tags.items(): + if seed_tag in index and tag in index: + res[(index[seed_tag], chain)] = index[tag] + return res + + def unroll_quartic_vertices(self, diaglist=None): + """Link every quartic vertex contribution to the cubic diagram it + merges with. + + Each colour structure of an unrollable quartic vertex (see + get_unrollable_quartic_vertices) splits its four legs into two pairs, + which is exactly what two cubic vertices joined by an internal line + do. Replacing every quartic vertex of a diagram by that pair of cubic + vertices therefore turns it into a diagram already present in the + amplitude, and the two carry the same colour factor up to a rational + coefficient. Their amplitudes can hence be summed before the colour + algebra is applied. + + Returns {(diagram_index, colour_chain): (target_index, coeff)} where + colour_chain follows the convention of color_amp.ColorBasis.colorize + (one colour structure index per vertex) and the amplitude of the + source is coeff times a contribution to the target diagram. Diagrams + without any quartic vertex are left out. + + To line the result up with HelasAmplitude, call this on the + *reconstructed* amplitude, HelasMatrixElement.get_base_amplitude(), + and not on the amplitude the matrix element was built from. That is + the one helas_objects itself colorizes, and the vertex order of the + reconstructed diagrams -- hence the position of a colour structure + index inside the chain -- can differ from the generated ones. The two + happen to agree up to six gluons and start to differ at seven. + """ + + model = self.get('process').get('model') + if diaglist is None: + diaglist = self.get('diagrams') + + unrollable = get_unrollable_quartic_vertices(model) + if not unrollable: + return {} + + quartic_positions = [] + for diag in diaglist: + quartic_positions.append([i for i, vx in + enumerate(diag.get('vertices')) + if vx.get('id') in unrollable]) + if not any(quartic_positions): + return {} + + ninitial = self.get_ninitial() + # Diagrams free of quartic vertices are the possible merge targets + target_index = {} + for i, diag in enumerate(diaglist): + if quartic_positions[i]: + continue + target_index[str(UnrollDiagramTag(diag, model, ninitial))] = i + + res = {} + for i, diag in enumerate(diaglist): + positions = quartic_positions[i] + if not positions: + continue + # only the colour structures which carry a coupling contribute, + # matching what color_amp.ColorBasis.colorize keeps + allowed = [sorted(misc.make_unique( + [key[0] for key in model.get_interaction( + diag.get('vertices')[p].get('id')).get('couplings')])) + for p in positions] + for keys in itertools.product(*allowed): + choice = dict(zip(positions, keys)) + unrolled = self.unrolled_diagram(diag, choice, unrollable) + try: + target = target_index[ + str(UnrollDiagramTag(unrolled, model, ninitial))] + except KeyError: + # No cubic partner: the merge is not available, which can + # happen when the partner was removed by a diagram filter + # or by a forbidden s-channel. + continue + chain = tuple(choice.get(p, 0) + for p in range(len(diag.get('vertices')))) + source_sig = diagram_colour_signature(diag, model, chain, + unrollable) + target_chain = (0,) * len(diaglist[target].get('vertices')) + target_sig = diagram_colour_signature(diaglist[target], model, + target_chain, unrollable) + if source_sig is None or target_sig is None or \ + source_sig[0] != target_sig[0]: + continue + res[(i, chain)] = (target, source_sig[1] / target_sig[1]) + + return res + + def unrolled_diagram(self, diagram, choice, unrollable): + """Return a copy of diagram where the quartic vertices listed in + choice (position -> colour structure index) have been replaced by the + two cubic vertices that colour structure factorises into.""" + + model = self.get('process').get('model') + vertices = base_objects.VertexList() + for i, vertex in enumerate(diagram.get('vertices')): + if i not in choice: + vertices.append(vertex) + continue + cubic_id, pairings = unrollable[vertex.get('id')] + vertices.extend(split_quartic_vertex( + vertex, i == len(diagram.get('vertices')) - 1, + pairings[choice[i]], cubic_id, model)) + + return base_objects.Diagram({'vertices': vertices}) + def apply_4gluon_specials(self, diag_list): - res = diag_list.__class__() + res = diag_list.__class__() for diag in diag_list: keep = True for vertex in diag.get('vertices'): @@ -988,9 +1668,14 @@ def copy_leglist(self, legs): [ copy.copy(leg) for leg in legs ]) def reduce_leglist(self, curr_leglist, max_multi_to1, ref_dict_to0, - is_decay_proc = False, coupling_orders = None): + is_decay_proc = False, coupling_orders = None, + cubic_legs = frozenset()): """Recursive function to reduce N LegList to N-1 For algorithm, see doc for generate_diagrams. + + cubic_legs holds the numbers of the legs of curr_leglist which were + produced by a vertex the seed rule keeps apart, and is only ever + non-empty when that rule is active. """ # Result variable which is a list of lists of vertices @@ -1023,6 +1708,8 @@ def reduce_leglist(self, curr_leglist, max_multi_to1, ref_dict_to0, vertex_id in vertex_ids] # Check for coupling orders. If orders < 0, skip vertex for final_vertex in final_vertices: + if self.joins_two_cubics(final_vertex, cubic_legs): + continue if self.reduce_orders(coupling_orders, model, [final_vertex.get('id')]) != False: res.append([final_vertex]) @@ -1060,13 +1747,24 @@ def reduce_leglist(self, curr_leglist, max_multi_to1, ref_dict_to0, # Some coupling order < 0 continue + # Seed rule: drop the combinations putting two of the cubic + # vertices to be kept apart on the same line + if any(self.shares_line_with_cubic(vertex.get('id'), + vertex.get('legs')[:-1], + cubic_legs) + for vertex in leg_vertex_tuple[1]): + continue + # This is where recursion happens # First, reduce again the leg part reduced_diagram = self.reduce_leglist(leg_vertex_tuple[0], max_multi_to1, ref_dict_to0, is_decay_proc, - new_coupling_orders) + new_coupling_orders, + self.mark_cubic_legs(\ + cubic_legs, + leg_vertex_tuple[1])) # If there is a reduced diagram if reduced_diagram: vertex_list_list = [list(leg_vertex_tuple[1])] @@ -1076,6 +1774,49 @@ def reduce_leglist(self, curr_leglist, max_multi_to1, ref_dict_to0, return res + def shares_line_with_cubic(self, vertex_id, incoming, cubic_legs): + """True when vertex_id is one of the cubic vertices the seed rule + keeps apart and one of the lines coming in was produced by another + one of them.""" + + if not self.seed_forbidden_cubic_ids or \ + vertex_id not in self.seed_forbidden_cubic_ids: + return False + return any(leg.get('number') in cubic_legs for leg in incoming) + + def joins_two_cubics(self, vertex, cubic_legs): + """True when the vertex closing the diagram puts two of the cubic + vertices the seed rule keeps apart on the same line. + + The closing vertex is either a real n->0 interaction, whose incoming + lines are simply the legs left over, or the identity vertex, which + states that its two legs are the two ends of one and the same line -- + and that line joins the two vertices which produced them.""" + + if not self.seed_forbidden_cubic_ids: + return False + legs = vertex.get('legs') + if vertex.get('id'): + return self.shares_line_with_cubic(vertex.get('id'), legs, + cubic_legs) + return all(leg.get('number') in cubic_legs for leg in legs) + + def mark_cubic_legs(self, cubic_legs, vertices): + """Update the set of leg numbers standing for a line produced by one + of the cubic vertices the seed rule keeps apart. + + A number is dropped as soon as the line is consumed, since a vertex + reuses the smallest number coming in for the leg it produces.""" + + if not self.seed_forbidden_cubic_ids: + return cubic_legs + consumed = frozenset(leg.get('number') for vertex in vertices + for leg in vertex.get('legs')[:-1]) + produced = frozenset(vertex.get('legs')[-1].get('number') + for vertex in vertices if vertex.get('id') in \ + self.seed_forbidden_cubic_ids) + return (cubic_legs - consumed) | produced + def reduce_orders(self, coupling_orders, model, vertex_id_list): """Return False if the coupling orders for any coupling is < 0, otherwise return the new coupling orders with the vertex diff --git a/madgraph/core/helas_objects.py b/madgraph/core/helas_objects.py index 46d3b015e..eaf4e649b 100755 --- a/madgraph/core/helas_objects.py +++ b/madgraph/core/helas_objects.py @@ -3997,6 +3997,15 @@ def default_setup(self): self._flavor_populated = False self._flavor_allow_trimming = False self._flavor_trimmed = False + # Cache for get_quartic_amplitude_merges(), needed both by the helas + # calls and by the colour amplitudes. Runtime only, like the above. + self.quartic_amplitude_merges = None + # Cache for get_quartic_current_sums(), same reason + self.quartic_current_sums = None + # Slots the current sums were given by reuse_outdated_wavefunctions + self.quartic_sum_me_ids = None + # Cache for get_amplitude_slots(), the recycled AMP array + self.amplitude_slots = None def filter(self, name, value): """Filter for valid diagram property values.""" @@ -4415,26 +4424,67 @@ def reuse_outdated_wavefunctions(self, helas_diagrams): for diag in helas_diagrams: for wf in diag['wavefunctions']: wf.set('me_id',wf.get('number')) + self.quartic_sum_me_ids = None # a fresh slot each, at the end return helas_diagrams + # A current sum is a line of its own, written as soon as the later of + # the two currents it reads is made, and read by the amplitudes it was + # built for. Giving it a key here lets it take a slot from the same + # pool as everything else, rather than one of its own at the end -- + # and lets the cubic current die at the sum rather than at the + # amplitude, since the amplitude no longer reads it. + sums, uses, folded = self.get_quartic_current_sums() + read_after = {} + for isum, (cubic, quartic, coeff) in enumerate(sums): + read_after.setdefault(max(cubic.get('number'), + quartic.get('number')), []).append(isum) + # keys for the sums, above every wavefunction number + offset = max([wf.get('number') for diag in helas_diagrams + for wf in diag['wavefunctions']] or [0]) + sum_key = lambda isum: offset + 1 + isum + # First compute the first/last appearance of each wavefunctions # first takes the line number and return the id of the created wf # last_lign takes the id of the wf and return the line number last_lign={} first={} pos=0 + written = set() + allocated = set() for diag in helas_diagrams: for wf in diag['wavefunctions']: pos+=1 for wfin in wf.get('mothers'): last_lign[wfin.get('number')] = pos assert wfin.get('number') in list(first.values()) + # the same wavefunction can be listed by more than one + # diagram; it is written twice with the same value, so it owns + # one slot from its first appearance to its last use, and + # handing it a second one here would leak the first + if wf.get('number') in allocated: + continue + allocated.add(wf.get('number')) first[pos] = wf.get('number') + for isum in read_after.get(wf.get('number'), []): + if isum in written: + continue + written.add(isum) + pos+=1 + cubic, quartic, coeff = sums[isum] + last_lign[cubic.get('number')] = pos + last_lign[quartic.get('number')] = pos + first[pos] = sum_key(isum) for amp in diag['amplitudes']: pos+=1 + substitution = uses.get(amp.get('number'), {}) for wfin in amp.get('mothers'): - last_lign[wfin.get('number')] = pos - + isum = substitution.get(wfin.get('number')) + if isum is None: + last_lign[wfin.get('number')] = pos + else: + # this amplitude reads the sum, not the cubic current + last_lign[sum_key(isum)] = pos + # last takes the line number and return the last appearing wf at #that particular line last=collections.defaultdict(list) @@ -4463,7 +4513,9 @@ def reuse_outdated_wavefunctions(self, helas_diagrams): for diag in helas_diagrams: for wf in diag['wavefunctions']: wf.set('me_id', replace[wf.get('number')]) - + self.quartic_sum_me_ids = [replace[sum_key(isum)] + for isum in range(len(sums))] + return helas_diagrams def restore_original_wavefunctions(self): @@ -4476,7 +4528,8 @@ def restore_original_wavefunctions(self): for diag in helas_diagrams: for wf in diag['wavefunctions']: wf.set('me_id',wf.get('number')) - + self.quartic_sum_me_ids = None # a fresh slot each, at the end + return helas_diagrams @@ -5314,12 +5367,15 @@ def get_num_configs(self): def get_number_of_wavefunctions(self): """Gives the total number of wavefunctions for this ME""" - out = max([wf.get('me_id') for wfs in self.get('diagrams') + # a current sum can hold the highest slot of all + extra = self.get_quartic_sum_me_ids() + + out = max([wf.get('me_id') for wfs in self.get('diagrams') for wf in wfs.get('wavefunctions')]) - if out: - return out - return sum([ len(d.get('wavefunctions')) for d in \ - self.get('diagrams')]) + if out: + return max([out] + extra) + return max([sum([ len(d.get('wavefunctions')) for d in \ + self.get('diagrams')])] + extra) def get_all_wavefunctions(self): """Gives a list of all wavefunctions for this ME""" @@ -6107,13 +6163,367 @@ def generate_color_amplitudes(self, color_basis, diagrams): return col_amp_list - def get_color_amplitudes(self): + + def get_color_amplitudes(self, merge_quartic_amplitudes=True): """Return a list of (coefficient, amplitude number) lists, corresponding to the JAMPs for this matrix element. The coefficients are given in the format (fermion factor, color - coeff (frac), imaginary, Nc power).""" - - return self.generate_color_amplitudes(self['color_basis'],self['diagrams']) + coeff (frac), imaginary, Nc power). + + merge_quartic_amplitudes says whether the caller also writes out the + sums of get_quartic_amplitude_merges. Only the Fortran writer does. + A backend which does not has to leave those amplitudes in the JAMPs, + where their own colour coefficients give the identical result -- the + two carry the same colour factor, which is the whole premise.""" + + col_amps = self.generate_color_amplitudes(self['color_basis'], + self['diagrams']) + # never computed at all: their current was summed instead + dropped = set(self.get_quartic_current_sums()[2]) + if merge_quartic_amplitudes: + # summed into their partner by GET_AMP, so they must not enter the + # JAMPs a second time + dropped |= set(self.get_quartic_amplitude_merges()) + if not dropped: + return col_amps + return [[entry for entry in col_amp if entry[1] not in dropped] + for col_amp in col_amps] + + def get_quartic_amplitude_merges(self): + """Return {amplitude number: (target number, coefficient)} for the + quartic gluon contributions which can be summed into another + amplitude. + + Each colour structure of a four gluon vertex carries the same colour + factor as the diagram obtained by splitting that vertex into two cubic + ones, see diagram_generation.Amplitude.unroll_quartic_vertices, so the + two amplitudes may be added before the colour algebra is applied. + Doing so here keeps the JAMPs -- and hence the work of the JAMP + optimiser -- proportional to the number of cubic diagrams rather than + to the number of amplitudes. + + The result is cached, since it is needed both when writing the helas + calls and when building the colour amplitudes. + """ + + if self.quartic_amplitude_merges is None: + self.quartic_amplitude_merges = self.compute_quartic_amplitude_merges() + return self.quartic_amplitude_merges + + def compute_quartic_amplitude_merges(self): + """Work out the amplitude merges, see get_quartic_amplitude_merges.""" + + if not madgraph.merge_quartic_vertices: + return {} + # colorize() is applied to the reconstructed amplitude, so this is the + # one whose colour chains line up with our 'color_indices' + base = self.get('base_amplitude') + links = base.unroll_quartic_vertices() + if not links: + return {} + + numbers = {} + for diagram in self.get('diagrams'): + for amplitude in diagram.get('amplitudes'): + numbers[(diagram.get('number') - 1, + tuple(amplitude.get('color_indices')))] = \ + amplitude.get('number') + + res = {} + for source, (target, coeff) in links.items(): + target_chain = (0,) * \ + len(base.get('diagrams')[target].get('vertices')) + try: + res[numbers[source]] = (numbers[(target, target_chain)], coeff) + except KeyError: + # no amplitude for one of the two: leave them alone + continue + + # a merged amplitude must never itself be merged away, otherwise the + # contributions it absorbed would be lost + targets = set(target for target, _ in res.values()) + return dict((source, value) for source, value in res.items() + if source not in targets) + + def get_quartic_current_sums(self): + """Return the current sums which take a quartic amplitude away. + + Where a quartic current and the cubic current carrying the same colour + factor feed the same vertex, the two amplitudes they give differ by + that one line and by nothing else. Summing the two currents once and + calling the amplitude on the sum therefore gets both contributions out + of a single call: + + TMP = W1 + c*W4 + CALL VVV1_0(..., TMP, AMP(t)) + + instead of one call for AMP(t) and one for the quartic AMP(s) which is + then added to it. The sum is one addition and is shared by every + amplitude which uses it, so it replaces as many calls as it has users. + + Only the amplitudes are treated. A sum sitting deeper would have to be + carried through every current above it, and those currents are shared + with diagrams which must not get the extra term -- see + docs/gluon-quartic-plan.md. + + Returns (sums, uses, folded): + sums [(cubic wavefunction, quartic wavefunction, coefficient)] + uses {amplitude number: {cubic wavefunction number: sums index}} + folded set of amplitude numbers the sums make unnecessary + """ + + if self.quartic_current_sums is None: + self.quartic_current_sums = self.compute_quartic_current_sums() + return self.quartic_current_sums + + def get_quartic_sum_me_ids(self): + """Wavefunction slot of each current sum. + + reuse_outdated_wavefunctions hands them out of the same pool as the + wavefunctions themselves, so a sum lands in whatever slot happens to + be free. When the wavefunctions were not recycled they get a fresh + slot each, at the end.""" + + sums = self.get_quartic_current_sums()[0] + if not sums: + return [] + if self.quartic_sum_me_ids is not None: + return self.quartic_sum_me_ids + used = [wf.get('me_id') or wf.get('number') + for diagram in self.get('diagrams') + for wf in diagram.get('wavefunctions')] + base = max(used or [0]) + return [base + 1 + isum for isum in range(len(sums))] + + def get_amplitude_slots(self): + """Recycle the AMP array the way reuse_outdated_wavefunctions recycles + the wavefunctions. + + Returns (slots, nslots, folds_at) where slots maps an amplitude number + onto its entry in AMP, nslots is how many entries that needs, and + folds_at maps a position in the emission order onto the merges to + write out just after it. + + An amplitude read by the JAMPs -- or by AMP2, which reads the same + ones, being the targets -- has to stay put until the end. A merge + source does not: once `AMP(t) = AMP(t) + AMP(s)` has run, its entry is + free. Writing each merge as soon as both of its amplitudes exist, + rather than all of them at the end, is what makes those entries worth + reclaiming. It only reclaims the gap between the two, so this is far + from the (2n-5)!! floor -- closing that would want the amplitudes + emitted in a different order, which is the wavefunction slot trade one + level down. + """ + + if self.amplitude_slots is not None: + return self.amplitude_slots + + folded = set(self.get_quartic_current_sums()[2]) + order = [amplitude.get('number') + for diagram in self.get('diagrams') + for amplitude in diagram.get('amplitudes') + if amplitude.get('number') not in folded] + position = dict((number, i) for i, number in enumerate(order)) + + # each merge is written as soon as both of its amplitudes are there + folds_at = {} + dies_at = {} + for source, (target, coeff) in \ + sorted(self.get_quartic_amplitude_merges().items()): + if source in folded or source not in position \ + or target not in position: + continue + at = max(position[source], position[target]) + folds_at.setdefault(at, []).append((target, source, coeff)) + dies_at.setdefault(at, []).append(source) + + slots, free, nslots = {}, [], 0 + for i, number in enumerate(order): + if free: + slots[number] = free.pop() + else: + nslots += 1 + slots[number] = nslots + # whatever this position's merges consume is free again after them + for source in dies_at.get(i, []): + free.append(slots[source]) + + self.amplitude_slots = (slots, nslots, folds_at) + return self.amplitude_slots + + def compute_quartic_current_sums(self): + """Work out the current sums, see get_quartic_current_sums.""" + + merges = self.get_quartic_amplitude_merges() + if not merges or madgraph.merge_quartic_vertices == 'slots': + # 'slots' orders the diagrams for the smallest wavefunction store + # instead, which puts a target ahead of the quartic currents + # feeding it, so no sum can be built + return [], {}, set() + + model = self.get('processes')[0].get('model') + unrollable = diagram_generation.get_unrollable_quartic_vertices(model) + cubic_ids = diagram_generation.get_unrollable_cubic_ids(model) + + amplitudes = {} + available = {} + written = set() + for diagram in self.get('diagrams'): + for wavefunction in diagram.get('wavefunctions'): + written.add(wavefunction.get('number')) + for amplitude in diagram.get('amplitudes'): + amplitudes[amplitude.get('number')] = amplitude + # the currents which exist by the time it is written out + available[amplitude.get('number')] = frozenset(written) + + # For each target, the substitutions each of its merges asks for + candidates = {} + for source, (target, coeff) in merges.items(): + pairs = self.match_quartic_mothers(amplitudes[source], + amplitudes[target], + unrollable, cubic_ids) + if not pairs: + continue + # the sum is written with sumw_1 / subw_1, which carry no weight + if abs(coeff) != 1: + continue + # the quartic current has to be there when the target is written + if any(quartic.get('number') not in available[target] + for cubic, quartic in pairs): + continue + key = frozenset((cubic.get('number'), quartic.get('number')) + for cubic, quartic in pairs) + candidates.setdefault(target, {})[key] = (source, coeff, pairs) + + sums, uses, folded = [], {}, set() + index = {} + for target in sorted(candidates): + entries = candidates[target] + singles = dict((next(iter(key)), value) + for key, value in entries.items() if len(key) == 1) + + # Substituting several mothers at once also produces the amplitude + # with all of them substituted, so every subset has to be a merge + # into this same target, with the product of the coefficients. + # Keep the substitutions which pass, drop the ones which do not. + chosen = [] + for pair in sorted(singles): + trial = chosen + [pair] + if all(self.subset_is_merged(entries, singles, combination) + for size in range(2, len(trial) + 1) + for combination in itertools.combinations(trial, size)): + chosen = trial + if not chosen: + continue + + for size in range(1, len(chosen) + 1): + for combination in itertools.combinations(chosen, size): + folded.add(entries[frozenset(combination)][0]) + for pair in chosen: + source, coeff, pairs = singles[pair] + cubic, quartic = pairs[0] + if (pair, coeff) not in index: + index[(pair, coeff)] = len(sums) + sums.append((cubic, quartic, coeff)) + uses.setdefault(target, {})[cubic.get('number')] = \ + index[(pair, coeff)] + + return sums, uses, folded + + @staticmethod + def subset_is_merged(entries, singles, combination): + """Is the amplitude with all of combination substituted a merge into + the same target, weighing the product of the single coefficients?""" + + entry = entries.get(frozenset(combination)) + if entry is None: + return False + weight = 1 + for pair in combination: + weight *= singles[pair][1] + return entry[1] == weight + + @staticmethod + def match_quartic_mothers(source, target, unrollable, cubic_ids): + """Pair the mothers up where the source amplitude carries the quartic + current and the target the cubic one taking the same four lines. + + Returns [] unless the two are otherwise one and the same vertex, which + is what makes the substitution a plain swap of one argument.""" + + if source.get('interaction_id') != target.get('interaction_id') or \ + source.get('color_key') != target.get('color_key'): + return [] + + source_mothers = dict((mother.get('number'), mother) + for mother in source.get('mothers')) + target_mothers = dict((mother.get('number'), mother) + for mother in target.get('mothers')) + only_source = [source_mothers[number] for number in source_mothers + if number not in target_mothers] + only_target = [target_mothers[number] for number in target_mothers + if number not in source_mothers] + if not only_source or len(only_source) != len(only_target): + return [] + + res = [] + for quartic in only_source: + hit = [cubic for cubic in only_target + if HelasMatrixElement.is_unrolled_pair(quartic, cubic, + unrollable, + cubic_ids)] + if len(hit) != 1: + return [] + res.append((hit[0], quartic)) + if len(set(cubic.get('number') for cubic, quartic in res)) != len(res): + return [] + return res + + @staticmethod + def is_unrolled_pair(quartic, cubic, unrollable, cubic_ids): + """True when the cubic current is the pair of vertices the quartic one + factorises into: same four lines coming in, same line going out, and + the colour structure the quartic carries is the one which separates + the two lines the inner cubic vertex joins. + + That last condition is the one which is easy to forget. A quartic + vertex makes one current per colour structure, and all of them take + the same lines and make the same line, so the lines alone cannot tell + them apart -- only the pairing can, and it has to be read against + sorted_mothers, the order ALOHA receives the legs in.""" + + if quartic.get('interaction_id') not in unrollable or \ + cubic.get('interaction_id') not in cubic_ids or \ + quartic.get('number_external') != cubic.get('number_external'): + return False + + pairings = unrollable[quartic.get('interaction_id')][1] + if quartic.get('color_key') >= len(pairings): + return False + pairing = pairings[quartic.get('color_key')] + mothers = HelasMatrixElement.sorted_mothers(quartic) + outgoing = quartic.find_outgoing_number() - 1 + slots = [i for i in range(len(mothers) + 1) if i != outgoing] + + wanted = sorted(mother.get('number') + for mother in quartic.get('mothers')) + for inner in cubic.get('mothers'): + if inner.get('interaction_id') not in cubic_ids or \ + len(inner.get('mothers')) != 2: + continue + if sorted([mother.get('number') + for mother in inner.get('mothers')] + + [other.get('number') for other in cubic.get('mothers') + if other is not inner]) != wanted: + continue + joined = set(mother.get('number') + for mother in inner.get('mothers')) + separated = frozenset(slots[i] for i, mother in enumerate(mothers) + if mother.get('number') in joined) + if frozenset(pairing[0]) == separated or \ + frozenset(pairing[1]) == separated: + return True + return False def sort_split_orders(self, split_orders): """ Sort the 'split_orders' list given in argument so that the orders of diff --git a/madgraph/fks/fks_base.py b/madgraph/fks/fks_base.py index 3b6f941d3..a16566d0f 100755 --- a/madgraph/fks/fks_base.py +++ b/madgraph/fks/fks_base.py @@ -129,6 +129,18 @@ def __init__(self, procdef=None, options={}): legs (stored in pdgs, so that they need to be generated only once and then reicycled """ + # The four gluon merging is validated at tree level only, and an NLO + # generation is left alone whatever the option says. fks reads the + # vertex decomposition of its real diagrams -- link_rb_configs finds + # the vertex splitting ij into i and j and takes it out -- and the + # unrolling re-roots them, which can put that pair in the closing + # vertex, where there is nothing to take out. + with misc.TMP_variable(madgraph, 'merge_quartic_vertices', False): + self.generate_all(procdef, options) + + def generate_all(self, procdef=None, options={}): + """Generates the born amplitudes, the born processes and the reals, + see __init__ which is the only caller.""" if 'nlo_mixed_expansion' in options: self['nlo_mixed_expansion'] = options['nlo_mixed_expansion'] diff --git a/madgraph/fks/fks_common.py b/madgraph/fks/fks_common.py index 7d055c6a0..c07305adb 100755 --- a/madgraph/fks/fks_common.py +++ b/madgraph/fks/fks_common.py @@ -196,14 +196,19 @@ def link_rb_configs(born_amp, real_amp, i, j, ij): for d in born_confs] - real_tags = [FKSDiagramTag(d['diagram'], - real_amp.get('process').get('model')) \ - for d in good_diags ] + # Dropping a duplicated tag has to drop its diagram with it: the two are + # walked in lockstep below, real_tags.index giving the position popped out + # of good_diags, so a tag left without its diagram misaligns everything + # after it. Which of two equal tags survives is decided by the order they + # come in, so keeping them apart also made the result order dependent. real_tags = [] + kept_diags = [] for d in good_diags: tag = FKSDiagramTag(d['diagram'], real_amp.get('process').get('model')) if not tag in real_tags: real_tags.append(tag) + kept_diags.append(d) + good_diags = kept_diags # and compare them if len(born_tags) != len(real_tags): diff --git a/madgraph/interface/madgraph_interface.py b/madgraph/interface/madgraph_interface.py index d7e1ef7fd..b1c8eee24 100755 --- a/madgraph/interface/madgraph_interface.py +++ b/madgraph/interface/madgraph_interface.py @@ -3155,7 +3155,8 @@ class MadGraphCmd(HelpToCmd, CheckValidForCmd, CompleteForCmd, CmdExtended): 'max_t_for_channel', 'zerowidth_tchannel', 'default_unset_couplings', - 'nlo_mixed_expansion' + 'nlo_mixed_expansion', + 'merge_quartic_vertices' ] _valid_nlo_modes = ['all','real','virt','sqrvirt','tree','noborn','LOonly', 'only'] _valid_sqso_types = ['==','<=','=','>'] @@ -3238,7 +3239,8 @@ class MadGraphCmd(HelpToCmd, CheckValidForCmd, CompleteForCmd, CmdExtended): 'max_t_for_channel': 99, # means no restrictions 'zerowidth_tchannel': True, 'nlo_mixed_expansion':True, - 'apply_flavor_grouping': True + 'apply_flavor_grouping': True, + 'merge_quartic_vertices': False } options_madevent = {'automatic_html_opening':True, @@ -3361,7 +3363,18 @@ def do_add(self, line): existing amplitudes or merge two model """ - + + # The four gluon merging is wanted while the diagrams are generated, + # which is below the interface, so it travels on the module. Synced + # here rather than only in the setter, since the option can also + # arrive from mg5_configuration.txt. + madgraph.merge_quartic_vertices = \ + self.options.get('merge_quartic_vertices', False) + # an added process arrives in the generated order whatever an earlier + # output reordered, so the two have to be brought back together -- + # a value matching neither makes the next output redo all of them + self._quartic_order = 'mixed' + args = self.split_arg(line) @@ -5007,6 +5020,10 @@ def clean_process(self): # Reset Helas matrix elements self._curr_matrix_elements = helas_objects.HelasMultiProcess() self._generate_info = "" + # Reset the diagrams kept for an 'auto' merge_quartic_vertices, they + # describe the amplitudes just dropped + self._quartic_order = None + self._quartic_pristine = None # Reset polarization-citation marker (a new process definition starts) self._uses_polarization = False self._uses_density_matrix = False @@ -9146,6 +9163,49 @@ def set2_zerowidth_tchannel(self, args, log=True): self.check_set(args) self.options[args[0]] = banner_module.ConfigFile.format_variable(args[1], bool, args[0]) + def help_set2_merge_quartic_vertices(self): + logger.info("merge_quartic_vertices ",'$MG:color:GREEN') + logger.info(" > (default: False) [pure gluon amplitudes]") + logger.info(" > Sum each four gluon contribution into the cubic") + logger.info(" > amplitude carrying the same colour factor.") + logger.info(" > False : off") + logger.info(" > speed : fewest amplitude calls (best on cpu)") + logger.info(" > slots : smallest wavefunction store (for gpu, where") + logger.info(" > that store is per thread); costs the sums") + logger.info(" > auto : below 6 legs only the amplitude merges, which") + logger.info(" > are free; above, slots when the matrix elements") + logger.info(" > go to a gpu backend and speed otherwise") + + def set2_merge_quartic_vertices(self, args, log=True): + """Sum the four gluon contributions into the cubic amplitude carrying + the same colour factor. + + Read while the diagrams are generated, so it has to be set before + 'generate' -- 'speed' and 'slots' fix the diagram order there and an + output time option would come too late. 'auto' generates in the + 'speed' order and lets each output reorder, see + apply_quartic_diagram_order. + + Example: set merge_quartic_vertices auto + """ + args = ['merge_quartic_vertices'] + args + self.check_set(args) + value = args[1].lower() + if value in ('slots', 'speed', 'auto'): + pass + elif value in ('false', '0', 'none', 'off'): + value = False + elif value in ('true', '1', 'on'): + value = 'speed' + else: + raise self.InvalidCmd( + "merge_quartic_vertices takes False, speed, slots or auto," + " not '%s'" % args[1]) + self.options[args[0]] = value + madgraph.merge_quartic_vertices = value + if log: + logger.info('set merge_quartic_vertices to %s' % value) + def set2_store_rwgt_info(self,args, log=True): """Set whether the code should generate systematics information in the output LHE file at NLO Default is set to False. @@ -9545,6 +9605,67 @@ def do_open(self, line): launch_ext.open_file(file_path) + # Output formats whose matrix elements can run on a gpu, where the + # wavefunction store is per thread. See set2_merge_quartic_vertices. + _gpu_me_formats = ['mg7', 'mg7_v5', 'standalone_mg7'] + # Diagram order currently materialised in _curr_amps, and the diagrams as + # they came out of the generation. Both only used for 'auto'. + _quartic_order = None + _quartic_pristine = None + + def apply_quartic_diagram_order(self, options): + """Resolve an 'auto' merge_quartic_vertices against the backend which + is about to be handed the matrix elements. + + What the two modes disagree on is the diagram order, and that is fixed + while the diagrams are generated. 'auto' generates in the 'speed' + order and reorders here instead, which is only possible because + 'slots' is that same order reversed. + + The reordering has to start from the diagrams as generated, not from + the ones in hand: an export mutates what it is given, and reversing + mutated diagrams gives an equivalent but differently numbered result. + """ + + if self.options.get('merge_quartic_vertices') != 'auto': + return + + target = options['me_exporter'].get('name', self._export_format) + gpu = target in self._gpu_me_formats or \ + options['me_exporter'].get('exporter', options['exporter']) == 'gpu' + wanted = 'slots' if gpu else 'speed' + madgraph.merge_quartic_vertices = wanted + + # Keep the diagrams of the amplitudes the seed rule applied to -- the + # order of any other one is not ours to touch. Both have to be read + # before the first export: it is the last moment the diagrams are the + # generated ones, and it also drops the marks saying they came from a + # seed. Amplitudes added later are picked up on the way past. + if self._quartic_pristine is None: + self._quartic_pristine = [] + for amp in self._curr_amps[len(self._quartic_pristine):]: + self._quartic_pristine.append( + copy.deepcopy(amp.get('diagrams')) + if (getattr(amp, 'seed_forbidden_cubic_ids', None) and + getattr(amp, 'quartic_unroll_tags', None)) else None) + if all(pristine is None for pristine in self._quartic_pristine): + return + logger.info("merge_quartic_vertices: '%s' for the %s matrix elements" + % (wanted, target)) + # the generation leaves them in the 'speed' order + if wanted == (self._quartic_order or 'speed'): + return + for amp, pristine in zip(self._curr_amps, self._quartic_pristine): + if pristine is None: + continue + diagrams = copy.deepcopy(pristine) + if wanted == 'slots': + diagrams.reverse() + amp.set('diagrams', base_objects.DiagramList(diagrams)) + self._quartic_order = wanted + # anything cached was built in the other order + self._curr_matrix_elements = helas_objects.HelasMultiProcess() + def do_output(self, line): """Main commands: Initialize a new Template or reinitialize one""" @@ -9668,7 +9789,11 @@ def do_output(self, line): options['me_exporter']['name'] = me_exporter else: options['me_exporter'] = {} - + + # now that the backend getting the matrix elements is known, an 'auto' + # merge_quartic_vertices can be resolved -- before anything is built + self.apply_quartic_diagram_order(options) + # check if os.path.realpath(self._export_dir) == os.getcwd(): if len(args) == 0: diff --git a/madgraph/interface/master_interface.py b/madgraph/interface/master_interface.py index b18de5252..0626731a1 100755 --- a/madgraph/interface/master_interface.py +++ b/madgraph/interface/master_interface.py @@ -612,6 +612,9 @@ def help_set2_include_lepton_initiated_processes(self, *args, **opts): def help_set2_loop_color_flows(self, *args, **opts): return self.cmd.help_set2_loop_color_flows(self, *args, **opts) + def help_set2_merge_quartic_vertices(self, *args, **opts): + return self.cmd.help_set2_merge_quartic_vertices(self, *args, **opts) + def help_set2_loop_optimized_output(self, *args, **opts): return self.cmd.help_set2_loop_optimized_output(self, *args, **opts) diff --git a/madgraph/iolibs/export_cpp.py b/madgraph/iolibs/export_cpp.py index 354815b36..bd26638f1 100755 --- a/madgraph/iolibs/export_cpp.py +++ b/madgraph/iolibs/export_cpp.py @@ -967,7 +967,8 @@ def get_process_class_definitions(self, write=True): replace_dict['nprocesses'] = self.nprocesses - color_amplitudes = self.matrix_elements[0].get_color_amplitudes() + color_amplitudes = self.matrix_elements[0].get_color_amplitudes( + merge_quartic_amplitudes=False) # Number of color flows replace_dict['ncolor'] = len(color_amplitudes) @@ -1045,7 +1046,7 @@ def get_process_function_definitions(self, write=True): # Extract process class name (for the moment same as file name) replace_dict['process_class_name'] = self.process_name - color_amplitudes = [me.get_color_amplitudes() for me in \ + color_amplitudes = [me.get_color_amplitudes(merge_quartic_amplitudes=False) for me in \ self.matrix_elements] replace_dict['initProc_lines'] = \ @@ -2124,7 +2125,7 @@ def get_process_function_definitions(self, write=True): # Extract process class name (for the moment same as file name) replace_dict['process_class_name'] = self.process_name - color_amplitudes = [me.get_color_amplitudes() for me in \ + color_amplitudes = [me.get_color_amplitudes(merge_quartic_amplitudes=False) for me in \ self.matrix_elements] replace_dict['initProc_lines'] = \ diff --git a/madgraph/iolibs/export_python.py b/madgraph/iolibs/export_python.py index 2b1ee5924..7b438024f 100755 --- a/madgraph/iolibs/export_python.py +++ b/madgraph/iolibs/export_python.py @@ -243,7 +243,10 @@ def get_jamp_lines(self, matrix_element): res_list = [] - for i, coeff_list in enumerate(matrix_element.get_color_amplitudes()): + # this writer emits no amplitude sums, so the quartic contributions + # have to stay in the JAMPs with their own colour coefficients + for i, coeff_list in enumerate(matrix_element.get_color_amplitudes( + merge_quartic_amplitudes=False)): res = "jamp[%d] = " % i diff --git a/madgraph/iolibs/export_v4.py b/madgraph/iolibs/export_v4.py index 525478ac5..e07ae7464 100755 --- a/madgraph/iolibs/export_v4.py +++ b/madgraph/iolibs/export_v4.py @@ -2200,6 +2200,30 @@ def get_multi_channel_dictionary(diagrams, config_map): return config_to_diag_dict + @staticmethod + def get_amplitude_slot_map(matrix_element): + """{amplitude number: AMP entry} when the AMP array is recycled, else + None. See HelasMatrixElement.get_amplitude_slots -- the writer emits + each amplitude into its entry, so everything reading AMP afterwards + has to go through the same map.""" + + if not isinstance(matrix_element, helas_objects.HelasMatrixElement): + return None + if not matrix_element.get_quartic_amplitude_merges(): + return None + return matrix_element.get_amplitude_slots()[0] + + @classmethod + def map_color_amplitudes(cls, matrix_element, color_amplitudes): + """The colour amplitudes with the amplitude numbers replaced by the + AMP entries they were written into.""" + + slots = cls.get_amplitude_slot_map(matrix_element) + if slots is None: + return color_amplitudes + return [[(coeff, slots[number]) for coeff, number in col_amp] + for col_amp in color_amplitudes] + def get_amp2_lines(self, matrix_element, config_map = [], replace_dict=None): """Return the amp2(i) = sum(amp for diag(i))^2 lines""" @@ -2227,7 +2251,10 @@ def get_amp2_lines(self, matrix_element, config_map = [], replace_dict=None): line = "AMP2(%(num)d)=AMP2(%(num)d)+" % \ {"num": (config_to_diag_dict[config][0] + 1)} - amp = "+".join(["AMP(%(num)d)" % {"num": a.get('number')} for a in \ + slots = self.get_amplitude_slot_map(matrix_element) + amp = "+".join(["AMP(%(num)d)" % + {"num": slots[a.get('number')] if slots + else a.get('number')} for a in \ sum([diagrams[idiag].get('amplitudes') for \ idiag in config_to_diag_dict[config]], [])]) @@ -2323,7 +2350,8 @@ def get_JAMP_lines_split_order(self, col_amps, split_order_amps, error_msg="Malformed '%s' argument passed to the "+\ "get_JAMP_lines_split_order function: %s"%str(split_order_amps) if(isinstance(col_amps,helas_objects.HelasMatrixElement)): - color_amplitudes=col_amps.get_color_amplitudes() + color_amplitudes=self.map_color_amplitudes( + col_amps, col_amps.get_color_amplitudes()) elif(isinstance(col_amps,list)): if(col_amps and isinstance(col_amps[0],list)): color_amplitudes=col_amps @@ -2397,7 +2425,8 @@ def get_JAMP_lines(self, col_amps, JAMP_format="JAMP(%s)", AMP_format="AMP(%s)", # Let the user call get_JAMP_lines directly from a MatrixElement or from # the color amplitudes lists. if(isinstance(col_amps,helas_objects.HelasMatrixElement)): - color_amplitudes=col_amps.get_color_amplitudes() + color_amplitudes=self.map_color_amplitudes( + col_amps, col_amps.get_color_amplitudes()) elif(isinstance(col_amps,list)): if(col_amps and isinstance(col_amps[0],list)): color_amplitudes=col_amps @@ -4240,7 +4269,10 @@ def write_matrix_element_v4(self, writer, matrix_element, fortran_model, # Extract ngraphs ngraphs = matrix_element.get_number_of_amplitudes() - replace_dict['ngraphs'] = ngraphs + # NGRAPHS only dimensions AMP, and AMP is recycled + slots = self.get_amplitude_slot_map(matrix_element) + replace_dict['ngraphs'] = \ + matrix_element.get_amplitude_slots()[1] if slots else ngraphs # Extract nwavefuncs nwavefuncs = matrix_element.get_number_of_wavefunctions() @@ -4601,7 +4633,8 @@ def get_JAMP_lines(self, col_amps, JAMP_format="JAMP(%s)", AMP_format="AMP(%s)", error_msg="Malformed '%s' argument passed to the get_JAMP_lines" if(isinstance(col_amps,helas_objects.HelasMatrixElement)): - col_amps=col_amps.get_color_amplitudes() + col_amps=self.map_color_amplitudes( + col_amps, col_amps.get_color_amplitudes()) elif(isinstance(col_amps,list)): if(col_amps and isinstance(col_amps[0],list)): col_amps=col_amps @@ -5128,7 +5161,10 @@ def write_matrix_element_v4(self, writer, matrix_element, fortran_model,proc_id # Extract ngraphs ngraphs = matrix_element.get_number_of_amplitudes() - replace_dict['ngraphs'] = ngraphs + # NGRAPHS only dimensions AMP, and AMP is recycled + slots = self.get_amplitude_slot_map(matrix_element) + replace_dict['ngraphs'] = \ + matrix_element.get_amplitude_slots()[1] if slots else ngraphs # Extract nwavefuncs nwavefuncs = matrix_element.get_number_of_wavefunctions() @@ -6174,7 +6210,10 @@ def write_matrix_element_v4(self, writer, matrix_element, fortran_model, # Extract ngraphs ngraphs = matrix_element.get_number_of_amplitudes() - replace_dict['ngraphs'] = ngraphs + # NGRAPHS only dimensions AMP, and AMP is recycled + slots = self.get_amplitude_slot_map(matrix_element) + replace_dict['ngraphs'] = \ + matrix_element.get_amplitude_slots()[1] if slots else ngraphs # Extract ndiags ndiags = len(matrix_element.get('diagrams')) diff --git a/madgraph/iolibs/helas_call_writers.py b/madgraph/iolibs/helas_call_writers.py index 2910e6a66..06f813f71 100755 --- a/madgraph/iolibs/helas_call_writers.py +++ b/madgraph/iolibs/helas_call_writers.py @@ -231,19 +231,129 @@ def get_matrix_element_calls(self, matrix_element): me = matrix_element.get('diagrams') matrix_element.reuse_outdated_wavefunctions(me) + # A current sum has to be written out once both currents are there. + # Doing it as soon as the later of the two is made keeps them alive: + # a slot is only handed on after its last use, and the cubic one is + # still needed by the amplitude the sum is for. + sums, uses, folded = self.get_quartic_current_sums(matrix_element) + slots = matrix_element.get_quartic_sum_me_ids() + after = {} + for i, (cubic, quartic, coeff) in enumerate(sums): + after.setdefault(max(cubic.get('number'), quartic.get('number')), + []).append(i) + + # The AMP array is recycled the same way, where the writer knows how + # to. A merge is written as soon as both of its amplitudes are there, + # which is what frees the source's entry, see get_amplitude_slots. + amp_slots = self.get_amplitude_slots(matrix_element) + position = [0] + res = [] + written = set() for diagram in matrix_element.get('diagrams'): - - - res.extend([ self.get_wavefunction_call(wf) for \ - wf in diagram.get('wavefunctions') ]) + + + for wf in diagram.get('wavefunctions'): + res.append(self.get_wavefunction_call(wf)) + for i in after.get(wf.get('number'), []): + # a wavefunction number can be listed by more than one + # diagram, and the sum must only be written once + if i in written: + continue + written.add(i) + cubic, quartic, coeff = sums[i] + res.extend(self.get_current_sum_lines( + slots[i], cubic, quartic, coeff)) res.append("# Amplitude(s) for diagram number %d" % \ diagram.get('number')) for amplitude in diagram.get('amplitudes'): - res.append(self.get_amplitude_call(amplitude)) + if amplitude.get('number') in folded: + # summed into another amplitude through a current sum + continue + res.append(self.get_amplitude_call_on_slot( + amplitude, uses.get(amplitude.get('number')), slots, + amp_slots)) + if amp_slots is not None: + res.extend(self.get_amplitude_merge_lines_at( + amp_slots, position[0])) + position[0] += 1 + + if amp_slots is None: + res.extend(self.get_amplitude_merge_lines(matrix_element)) return res + def get_amplitude_slots(self, matrix_element): + """The recycled AMP array, or None to leave AMP indexed by amplitude + number. Only the Fortran writer has an AMP array to recycle.""" + + return None + + def get_amplitude_call_on_slot(self, amplitude, substitution, slots, + amp_slots): + """The amplitude call, written into its recycled AMP entry. + + The slot is swapped onto the amplitude's number and put straight back, + the same way get_amplitude_call_on_sums does it for the mothers: the + call is formatted from `out`, which is that number.""" + + if amp_slots is None: + return self.get_amplitude_call_on_sums(amplitude, substitution, + slots) + number = amplitude.get('number') + amplitude.set('number', amp_slots[0][number]) + try: + return self.get_amplitude_call_on_sums(amplitude, substitution, + slots) + finally: + amplitude.set('number', number) + + def get_amplitude_merge_lines_at(self, amp_slots, position): + """The merges due just after this amplitude. Fortran only.""" + + return [] + + def get_quartic_current_sums(self, matrix_element): + """The current sums to write out. Only the Fortran writer knows how to + emit one, see FortranUFOHelasCallWriter.""" + + return [], {}, set() + + def get_current_sum_lines(self, number, cubic, quartic, coeff): + """Lines building one current sum. Fortran only.""" + + raise NotImplementedError + + def get_amplitude_call_on_sums(self, amplitude, substitution, slots): + """The amplitude call, reading the current sums in place of the cubic + currents they were built from. + + The slot is swapped on the mother itself and put back straight away, + the same way get_loop_amplitude_helas_calls relabels its externals.""" + + if not substitution: + return self.get_amplitude_call(amplitude) + + original = [] + for mother in amplitude.get('mothers'): + index = substitution.get(mother.get('number')) + if index is None: + continue + original.append((mother, mother.get('me_id'))) + mother.set('me_id', slots[index]) + try: + return self.get_amplitude_call(amplitude) + finally: + for mother, me_id in original: + mother.set('me_id', me_id) + + def get_amplitude_merge_lines(self, matrix_element): + """Lines summing the quartic contributions into the amplitude which + carries the same colour factor. Only the Fortran writer implements + this, see FortranUFOHelasCallWriter.""" + + return [] + def get_wavefunction_calls(self, wavefunctions): """Return a list of strings, corresponding to the Helas calls for the matrix element""" @@ -1027,6 +1137,86 @@ class FortranUFOHelasCallWriter(UFOHelasCallWriter): mp_prefix = check_param_card.ParamCard.mp_prefix + def get_amplitude_merge_lines(self, matrix_element): + """Sum every quartic contribution into the amplitude carrying the + same colour factor. + + The two share a colour factor up to a rational coefficient, so adding + them here lets the JAMPs run over the cubic diagrams only, which is + what the JAMP optimiser then has to work with. The amplitudes summed + away are dropped from the colour amplitudes by + HelasMatrixElement.get_color_amplitudes.""" + + merges = matrix_element.get_quartic_amplitude_merges() + if not merges: + return [] + # these were never computed: their current was summed instead + folded = self.get_quartic_current_sums(matrix_element)[2] + + res = ['# Sum the quartic contributions into their cubic partner'] + for source in sorted(merges): + if source in folded: + continue + target, coeff = merges[source] + if coeff == 1: + res.append('AMP(%d) = AMP(%d) + AMP(%d)' % + (target, target, source)) + elif coeff == -1: + res.append('AMP(%d) = AMP(%d) - AMP(%d)' % + (target, target, source)) + else: + res.append('AMP(%d) = AMP(%d) + (%.15e)*AMP(%d)' % + (target, target, float(coeff), source)) + return res + + def get_amplitude_slots(self, matrix_element): + """The recycled AMP array, see + HelasMatrixElement.get_amplitude_slots.""" + + if not matrix_element.get_quartic_amplitude_merges(): + return None + return matrix_element.get_amplitude_slots() + + def get_amplitude_merge_lines_at(self, amp_slots, position): + """The merges whose two amplitudes are both there as of this + position, written into the slots they were given.""" + + slots, nslots, folds_at = amp_slots + res = [] + for target, source, coeff in folds_at.get(position, []): + args = (slots[target], slots[target], slots[source]) + if coeff == 1: + res.append('AMP(%d) = AMP(%d) + AMP(%d)' % args) + elif coeff == -1: + res.append('AMP(%d) = AMP(%d) - AMP(%d)' % args) + else: + res.append('AMP(%d) = AMP(%d) + (%.15e)*AMP(%d)' % + (slots[target], slots[target], float(coeff), + slots[source])) + return res + + def get_quartic_current_sums(self, matrix_element): + """The current sums, see HelasMatrixElement.get_quartic_current_sums""" + + return matrix_element.get_quartic_current_sums() + + def get_current_sum_lines(self, number, cubic, quartic, coeff): + """Sum the quartic current into the cubic one carrying the same colour + factor, so that the amplitude reading the sum gets both at once. + + Written as a call rather than as two assignments so that everything + reading these files -- the helicity recycling in particular, which + rebuilds the DAG from the calls alone -- sees an ordinary internal + wavefunction taking two mothers. sumw_1 and subw_1 live in + aloha_functions.f; the coefficient is always +-1, see + HelasMatrixElement.compute_quartic_current_sums.""" + + return ['CALL %s(%s,%s,%s)' % ( + 'SUMW_1' if coeff == 1 else 'SUBW_1', + self.format_helas_object('W(', '%d') % cubic.get('me_id'), + self.format_helas_object('W(', '%d') % quartic.get('me_id'), + self.format_helas_object('W(', '%d') % number)] + def __init__(self, argument={}, hel_sum = False, options={}): """Allow generating a HelasCallWriter from a Model.The hel_sum argument specifies if amplitude and wavefunctions must be stored specifying the diff --git a/madgraph/iolibs/template_files/madmatrix/cpp_hel_amps_h.inc b/madgraph/iolibs/template_files/madmatrix/cpp_hel_amps_h.inc index 8dde0fff4..bed053192 100644 --- a/madgraph/iolibs/template_files/madmatrix/cpp_hel_amps_h.inc +++ b/madgraph/iolibs/template_files/madmatrix/cpp_hel_amps_h.inc @@ -48,6 +48,39 @@ namespace mg5amcCpu : pvec(pvec_sv), w(reinterpret_cast(w_sv)), flv_index(flv) {} }; + // Sum two currents standing for the same off shell line: the four gluon + // current and the pair of three gluon vertices it factorises into carry the + // same colour factor, so the amplitude reading the sum gets both + // contributions from a single call. See + // HelasMatrixElement.get_quartic_current_sums. The two share their momentum, + // so only the wavefunction is added and the rest is taken from the first. + template + __device__ inline void + SUMW_1( const ALOHAOBJ& V2, const ALOHAOBJ& V3, ALOHAOBJ& V1 ) + { + const cxtype_sv* wV2 = W_ACCESS::kernelAccessConst( V2.w ); + const cxtype_sv* wV3 = W_ACCESS::kernelAccessConst( V3.w ); + cxtype_sv* wV1 = W_ACCESS::kernelAccess( V1.w ); + for( int i = 0; i < ALOHAOBJ::np4; i++ ) V1.pvec[i] = V2.pvec[i]; + for( int i = 0; i < ALOHAOBJ::nw6; i++ ) wV1[i] = wV2[i] + wV3[i]; + V1.flv_index = V2.flv_index; + return; + } + + // As SUMW_1, for the contributions which enter with a minus sign. + template + __device__ inline void + SUBW_1( const ALOHAOBJ& V2, const ALOHAOBJ& V3, ALOHAOBJ& V1 ) + { + const cxtype_sv* wV2 = W_ACCESS::kernelAccessConst( V2.w ); + const cxtype_sv* wV3 = W_ACCESS::kernelAccessConst( V3.w ); + cxtype_sv* wV1 = W_ACCESS::kernelAccess( V1.w ); + for( int i = 0; i < ALOHAOBJ::np4; i++ ) V1.pvec[i] = V2.pvec[i]; + for( int i = 0; i < ALOHAOBJ::nw6; i++ ) wV1[i] = wV2[i] - wV3[i]; + V1.flv_index = V2.flv_index; + return; + } + struct FLV_COUPLING_VIEW { const int* const partner1; diff --git a/madgraph/madevent/hel_recycle.py b/madgraph/madevent/hel_recycle.py index af0a1e7da..4dabf2556 100755 --- a/madgraph/madevent/hel_recycle.py +++ b/madgraph/madevent/hel_recycle.py @@ -497,8 +497,12 @@ def add_amp_index(self, matchobj): def add_indices(self, line): '''Add loop_var index to amp and output variable. Also update name of output variable.''' - # Doesnt work if the AMP arguments contain brackets - new_line = re.sub(r'\WAMP\(.*?\)', self.add_amp_index, line) + # Doesnt work if the AMP arguments contain brackets. + # The character in front is looked at rather than eaten, so that an + # AMP( opening the statement is indexed too -- which is what a line + # like "AMP(31) = AMP(31) + AMP(1)" needs. + new_line = re.sub(r'(?> iflavor ) & 0x1ULL ) {' % group_mask + # OM - the four gluon optimisation (set merge_quartic_vertices). A quartic + # current and the cubic current carrying the same colour factor are + # summed into a third one, which the amplitude reads instead, so that + # one call gets both contributions and the quartic amplitude is never + # computed. The sum is written as soon as the later of the two + # currents is made. Unlike Fortran there is no AMP array here, so the + # amplitudes which cannot be reached this way are simply left alone: + # their own colour coefficients put them in the right JAMPs, which is + # why get_color_amplitudes is asked not to drop them. + sums, sum_uses, sum_folded = matrix_element.get_quartic_current_sums() + sum_slots = matrix_element.get_quartic_sum_me_ids() + sum_written = set() + sum_after = {} + for isum, (cubic, quartic, coeff) in enumerate(sums): + sum_after.setdefault(max(cubic.get('number'), + quartic.get('number')), []).append(isum) + id_amp = 0 for diagram in matrix_element.get('diagrams'): ###print('DIAGRAM %3d: #wavefunctions=%3d, #diagrams=%3d' % @@ -2547,13 +2565,39 @@ def _guard_open(group_mask): res.append('}') else: res.append(call) + for isum in sum_after.get(wf.get('number'), []): + # a wavefunction number can be listed by more than one + # diagram here, and the sum must only be written once + if isum in sum_written: + continue + sum_written.add(isum) + cubic, quartic, coeff = sums[isum] + res.append('%s( aloha_obj[%d], aloha_obj[%d],' + ' aloha_obj[%d] );' + % ('SUMW_1' if coeff == 1 else 'SUBW_1', + cubic.get('me_id') - 1, + quartic.get('me_id') - 1, + sum_slots[isum] - 1)) if len(diagram.get('wavefunctions')) == 0 : res.append('// (none)') # AV res.append('\n // Amplitude(s) for diagram number %d' % diagram.get('number')) for amplitude in diagram.get('amplitudes'): id_amp +=1 + if amplitude.get('number') in sum_folded: + continue # summed into another amplitude as a current namp = amplitude.get('number') amplitude.set('number', 1) + # OM - read the current sum in place of the cubic current it + # was built from, the same way the Fortran writer does + sum_original = [] + for mother in amplitude.get('mothers'): + isum = sum_uses.get(namp, {}).get(mother.get('number')) + if isum is None: + continue + sum_original.append((mother, mother.get('me_id'))) + mother.set('me_id', sum_slots[isum]) amp_block = [ self.get_amplitude_call(amplitude) ] # AV new: avoid format_call + for mother, me_id in sum_original: + mother.set('me_id', me_id) if id_amp in diag_to_config: ###res.append("if( channelId == %i ) numerators_sv += cxabs2( amp_sv[0] );" % diag_to_config[id_amp]) # BUG #472 ###res.append("if( channelId == %i ) numerators_sv += cxabs2( amp_sv[0] );" % id_amp) # wrong fix for BUG #472 diff --git a/tests/unit_tests/core/test_color_amp.py b/tests/unit_tests/core/test_color_amp.py index 9775c450f..0926cae2c 100755 --- a/tests/unit_tests/core/test_color_amp.py +++ b/tests/unit_tests/core/test_color_amp.py @@ -20,6 +20,8 @@ import copy import fractions +import madgraph +import madgraph.various.misc as misc import madgraph.core.base_objects as base_objects import madgraph.core.diagram_generation as diagram_generation @@ -255,7 +257,11 @@ def test_colorize_uux_ggg(self): myamplitude.set('process', myprocess) - myamplitude.generate_diagrams() + # What is checked below is colorize, against diagrams picked by their + # position, so it wants the plain generation order -- the four gluon + # merging reorders them + with misc.TMP_variable(madgraph, 'merge_quartic_vertices', False): + myamplitude.generate_diagrams() my_col_basis = color_amp.ColorBasis() diff --git a/tests/unit_tests/core/test_diagram_generation.py b/tests/unit_tests/core/test_diagram_generation.py index 729c50ce3..954027b88 100755 --- a/tests/unit_tests/core/test_diagram_generation.py +++ b/tests/unit_tests/core/test_diagram_generation.py @@ -17,6 +17,7 @@ from __future__ import absolute_import import copy +import fractions import itertools import logging import math @@ -24,7 +25,10 @@ import tests.unit_tests as unittest +import madgraph +import madgraph.various.misc as misc import madgraph.core.base_objects as base_objects +import madgraph.core.color_amp as color_amp import madgraph.core.diagram_generation as diagram_generation import models.import_ufo as import_ufo from madgraph import MadGraph5Error, InvalidCmd @@ -3745,7 +3749,11 @@ def test_diagram_tag_gg_ggg(self): myproc = base_objects.Process({'legs':myleglist, 'model':self.base_model}) - myamplitude = diagram_generation.Amplitude(myproc) + # DiagramTag is what is checked here, against diagram numbers, so it + # wants the plain generation order -- the four gluon merging reorders + # them + with misc.TMP_variable(madgraph, 'merge_quartic_vertices', False): + myamplitude = diagram_generation.Amplitude(myproc) tags = [] permutations = [] @@ -3893,3 +3901,533 @@ def test_diagram_tag_to_diagram_uux_nglue(self): self.assertEqual(dtag, diagram_generation.DiagramTag(\ dtag.diagram_from_tag(self.base_model))) + +#=============================================================================== +# TestQuarticUnrolling +#=============================================================================== +class TestQuarticUnrolling(unittest.TestCase): + """Test the unrolling of quartic vertices into pairs of cubic ones""" + + def setUp(self): + self.base_model = import_ufo.import_model('sm') + + def make_amplitude(self, initial, final, orders=None): + myleglist = base_objects.LegList( + [base_objects.Leg({'id':pdg, 'state':False}) for pdg in initial] + + [base_objects.Leg({'id':pdg, 'state':True}) for pdg in final]) + mydict = {'legs':myleglist, 'model':self.base_model} + if orders: + mydict['orders'] = orders + return diagram_generation.Amplitude(base_objects.Process(mydict)) + + def colour_directions(self, amplitude): + """Return {(diagram, colour chain): (direction, norm)} where direction + is the colour vector of that contribution normalised to its first non + zero entry. Two contributions can be summed exactly when they share a + direction, the relative weight being the ratio of their norms. This is + computed from the colour algebra and is completely independent from + the unrolling being tested.""" + + basis = color_amp.ColorBasis() + basis.build(amplitude) + components = {} + for index, key in enumerate(sorted(basis.keys())): + for (diag, chain, coeff, imag, nc, loop_nc) in basis[key]: + components.setdefault((diag, chain), {})[index] = \ + (fractions.Fraction(coeff), imag, nc, loop_nc) + + res = {} + for piece, entries in components.items(): + first = min(entries) + norm, imag, nc, loop_nc = entries[first] + res[piece] = (tuple(sorted( + (i, c / norm, m ^ imag, n - nc, l - loop_nc) + for i, (c, m, n, l) in entries.items())), norm) + return res + + def check_process(self, initial, final, nlink, orders=None): + """Every link must reproduce the colour algebra, and every quartic + contribution must be linked to exactly one cubic diagram.""" + + amplitude = self.make_amplitude(initial, final, orders) + links = amplitude.unroll_quartic_vertices() + directions = self.colour_directions(amplitude) + diagrams = amplitude.get('diagrams') + + self.assertEqual(len(links), nlink) + + for (diag, chain), (target, coeff) in links.items(): + target_chain = (0,) * len(diagrams[target].get('vertices')) + source_dir, source_norm = directions[(diag, chain)] + target_dir, target_norm = directions[(target, target_chain)] + # same colour direction, and the coefficient is the relative weight + self.assertEqual(source_dir, target_dir) + self.assertEqual(coeff, source_norm / target_norm) + # the target is a genuine cubic diagram + self.assertEqual(amplitude.unrolled_diagram( + diagrams[target], {}, + diagram_generation.get_unrollable_quartic_vertices( + self.base_model)).get('vertices'), + diagrams[target].get('vertices')) + + # nothing is left behind: the contributions which are not linked are + # exactly the ones carried by a diagram without any quartic vertex + unrollable = diagram_generation.get_unrollable_quartic_vertices( + self.base_model) + cubic = [i for i, d in enumerate(diagrams) + if not any(v.get('id') in unrollable + for v in d.get('vertices'))] + unlinked = [piece for piece in directions if piece not in links] + self.assertEqual(sorted(piece[0] for piece in unlinked), sorted(cubic)) + + def test_unrollable_vertices_sm(self): + """The four gluon vertex is the only one factorising in the SM""" + + unrollable = diagram_generation.get_unrollable_quartic_vertices( + self.base_model) + self.assertEqual(len(unrollable), 1) + (cubic, pairings), = unrollable.values() + quartic, = unrollable.keys() + self.assertEqual([p.get_pdg_code() for p in + self.base_model.get_interaction(quartic).get('particles')], + [21, 21, 21, 21]) + self.assertEqual([p.get_pdg_code() for p in + self.base_model.get_interaction(cubic).get('particles')], + [21, 21, 21]) + # each colour structure splits the four legs into two pairs + self.assertEqual(pairings, [((0, 1), (2, 3)), + ((0, 2), (1, 3)), + ((0, 3), (1, 2))]) + + def test_unroll_gg_gg(self): + """g g > g g: the contact term merges into the s, t and u channel""" + + self.check_process([21, 21], [21, 21], 3) + + def test_unroll_gg_ggg(self): + """g g > g g g: 45 contributions collapse onto the 15 cubic diagrams""" + + self.check_process([21, 21], [21, 21, 21], 30) + + def test_unroll_gg_gggg(self): + """g g > g g g g: 510 contributions collapse onto 105 cubic diagrams""" + + self.check_process([21, 21], [21, 21, 21, 21], 405) + + def test_unroll_bbx_ggg(self): + """A quartic vertex sitting next to a colour triplet line""" + + self.check_process([5, -5], [21, 21, 21], 3) + + def test_unroll_gg_ttxg(self): + """The quartic vertex feeding an internal gluon of a t t~ pair""" + + self.check_process([21, 21], [6, -6, 21], 3) + + def test_links_match_helas_colour_indices(self): + """The links must be usable against HelasAmplitude. + + helas_objects colorizes the reconstructed base amplitude, whose + vertex order can differ from the generated diagrams, so the chains + only line up with HelasAmplitude.get('color_indices') when the + unrolling is run on get_base_amplitude().""" + + import madgraph.core.helas_objects as helas_objects + + amplitude = self.make_amplitude([21, 21], [21, 21, 21]) + matrix_element = helas_objects.HelasMatrixElement(amplitude) + base = matrix_element.get('base_amplitude') + + known = {} + for diagram in matrix_element.get('diagrams'): + for helas_amp in diagram.get('amplitudes'): + known[(diagram.get('number') - 1, + tuple(helas_amp.get('color_indices')))] = \ + helas_amp.get('number') + + links = base.unroll_quartic_vertices() + self.assertTrue(links) + for (diag, chain), (target, _) in links.items(): + target_chain = (0,) * len(base.get('diagrams')[target].get('vertices')) + self.assertIn((diag, chain), known) + self.assertIn((target, target_chain), known) + # every amplitude is either a merge target or folded into one + self.assertEqual(len(known) - len(links), + len(set(target for target, _ in links.values()))) + + def test_no_unrolling_without_quartic(self): + """Processes without a four gluon vertex give no link""" + + amplitude = self.make_amplitude([5, -5], [5, -5]) + self.assertEqual(amplitude.unroll_quartic_vertices(), {}) + +#=============================================================================== +# TestSeedRule +#=============================================================================== +class TestSeedRule(unittest.TestCase): + """Test the seed rule: no two cubic gluon vertices sharing a line. + + Unrolling a quartic vertex always yields two cubic vertices joined by the + line which replaced it, so the diagrams generation has to keep are exactly + those with no such pair to contract back.""" + + def setUp(self): + self.base_model = import_ufo.import_model('sm') + self.cubic_ids = diagram_generation.get_unrollable_cubic_ids( + self.base_model) + self.merge_quartic = madgraph.merge_quartic_vertices + + def tearDown(self): + madgraph.merge_quartic_vertices = self.merge_quartic + + def make_diagrams(self, initial, final, seed): + """The generated diagrams, either the full set or -- with the seed + rule on and the expansion stopped -- the seed it starts from.""" + + madgraph.merge_quartic_vertices = seed + myleglist = base_objects.LegList( + [base_objects.Leg({'id':pdg, 'state':False}) for pdg in initial] + + [base_objects.Leg({'id':pdg, 'state':True}) for pdg in final]) + process = base_objects.Process({'legs':myleglist, + 'model':self.base_model}) + if not seed: + return diagram_generation.Amplitude(process).get('diagrams') + + class SeedOnlyAmplitude(diagram_generation.Amplitude): + """Stops after the seed, so that it can be looked at""" + def expand_seed_diagrams(self, seed): + return seed + + return SeedOnlyAmplitude(process).get('diagrams') + + def cubic_adjacencies(self, diagram): + """Number of lines joining two cubic gluon vertices, counted without + any help from the generation. A line is recognised by its leg number + being live, a vertex reusing the smallest number coming in for the leg + it produces.""" + + vertices = diagram.get('vertices') + last = len(vertices) - 1 + live = {} + count = 0 + for i, vertex in enumerate(vertices): + legs = vertex.get('legs') + for leg in (legs if i == last else legs[:-1]): + producer = live.pop(leg.get('number'), None) + if producer is not None and \ + vertices[producer].get('id') in self.cubic_ids and \ + vertex.get('id') in self.cubic_ids: + count += 1 + if i != last: + live[legs[-1].get('number')] = i + self.assertFalse(live) + return count + + def check_process(self, initial, final, nfull, nseed): + """The generated seed has to be exactly the full generation filtered + on the rule -- same diagrams, not merely the same count.""" + + full = self.make_diagrams(initial, final, False) + seed = self.make_diagrams(initial, final, True) + self.assertEqual(len(full), nfull) + self.assertEqual(len(seed), nseed) + + def tags(diagrams): + return set(str(diagram_generation.UnrollDiagramTag( + diagram, self.base_model, len(initial))) + for diagram in diagrams) + + self.assertEqual(tags(seed), + tags([d for d in full + if not self.cubic_adjacencies(d)])) + self.assertFalse([d for d in seed if self.cubic_adjacencies(d)]) + + def test_seed_gg_gg(self): + """g g > g g: only the contact term has no cubic pair""" + + self.check_process([21, 21], [21, 21], 4, 1) + + def test_seed_gg_ggg(self): + """g g > g g g: the ten one quartic one cubic diagrams""" + + self.check_process([21, 21], [21, 21, 21], 25, 10) + + def test_seed_gg_gggg(self): + """g g > g g g g: 45 with a quartic in the middle, 10 with two""" + + self.check_process([21, 21], [21, 21, 21, 21], 220, 55) + + def test_seed_gg_ggggg(self): + """g g > g g g g g""" + + self.check_process([21, 21], [21, 21, 21, 21, 21], 2485, 385) + + def test_seed_gg_ttxgg(self): + """A cubic gluon vertex next to a quark line is not touched""" + + self.check_process([21, 21], [6, -6, 21, 21], 123, 84) + + def check_expansion(self, initial, final, nfull, nlink): + """Expanding the seed has to give the baseline diagram set back, and + the links it records have to be the ones the colour algebra gives.""" + + base = self.make_diagrams(initial, final, False) + madgraph.merge_quartic_vertices = True + myleglist = base_objects.LegList( + [base_objects.Leg({'id':pdg, 'state':False}) for pdg in initial] + + [base_objects.Leg({'id':pdg, 'state':True}) for pdg in final]) + amplitude = diagram_generation.Amplitude(base_objects.Process( + {'legs':myleglist, 'model':self.base_model})) + expanded = amplitude.get('diagrams') + + def tags(diagrams): + return set(str(diagram_generation.UnrollDiagramTag( + diagram, self.base_model, len(initial))) + for diagram in diagrams) + + # same diagrams, and no diagram reached twice + self.assertEqual(len(expanded), nfull) + self.assertEqual(len(base), nfull) + self.assertEqual(tags(expanded), tags(base)) + self.assertEqual(len(tags(expanded)), nfull) + + # the links recorded while expanding, and the independent ones + recorded = amplitude.get_quartic_unroll_links() + colour = amplitude.unroll_quartic_vertices() + self.assertEqual(len(recorded), nlink) + self.assertEqual(set(recorded), set(colour)) + for key, target in recorded.items(): + self.assertEqual(target, colour[key][0]) + + def test_expand_gg_gg(self): + """g g > g g: the contact term unrolls into s, t and u""" + + self.check_expansion([21, 21], [21, 21], 4, 3) + + def test_expand_gg_ggg(self): + """g g > g g g""" + + self.check_expansion([21, 21], [21, 21, 21], 25, 30) + + def test_expand_gg_gggg(self): + """g g > g g g g: 55 seeds reach all 220 diagrams""" + + self.check_expansion([21, 21], [21, 21, 21, 21], 220, 405) + + def test_expand_gg_ttxgg(self): + """Expansion leaves the quark lines alone""" + + self.check_expansion([21, 21], [6, -6, 21, 21], 123, 54) + + def test_current_sums_gg_ggg(self): + """g g > g g g: seven quartic amplitudes become a current sum""" + + self.check_current_sums([21, 21], [21, 21, 21], 7, 7) + + def test_current_sums_gg_gggg(self): + """g g > g g g g: thirty sums take sixty amplitude calls away""" + + self.check_current_sums([21, 21], [21, 21, 21, 21], 30, 60) + + def test_current_sums_inactive_by_default(self): + """No current sum unless madgraph.merge_quartic_vertices is set""" + + import madgraph.core.helas_objects as helas_objects + + madgraph.merge_quartic_vertices = False + myleglist = base_objects.LegList( + [base_objects.Leg({'id':21, 'state':False})] * 2 + + [base_objects.Leg({'id':21, 'state':True})] * 3) + matrix_element = helas_objects.HelasMatrixElement( + diagram_generation.Amplitude(base_objects.Process( + {'legs':myleglist, 'model':self.base_model}))) + self.assertEqual(matrix_element.get_quartic_current_sums(), + ([], {}, set())) + self.assertEqual(matrix_element.get_quartic_sum_me_ids(), []) + + def check_current_sums(self, initial, final, nsum, nfolded): + """A current sum has to stand for exactly the amplitude it takes away: + the same vertex, with the quartic current where the target has the + cubic one.""" + + import madgraph.core.helas_objects as helas_objects + + madgraph.merge_quartic_vertices = True + myleglist = base_objects.LegList( + [base_objects.Leg({'id':pdg, 'state':False}) for pdg in initial] + + [base_objects.Leg({'id':pdg, 'state':True}) for pdg in final]) + matrix_element = helas_objects.HelasMatrixElement( + diagram_generation.Amplitude(base_objects.Process( + {'legs':myleglist, 'model':self.base_model}))) + + sums, uses, folded = matrix_element.get_quartic_current_sums() + merges = matrix_element.get_quartic_amplitude_merges() + self.assertEqual(len(sums), nsum) + self.assertEqual(len(folded), nfolded) + # every sum gets a wavefunction slot, out of the same pool as the + # wavefunctions themselves + slots = matrix_element.get_quartic_sum_me_ids() + self.assertEqual(len(slots), nsum) + self.assertTrue(all(slot > 0 for slot in slots)) + + amplitudes = dict((amplitude.get('number'), amplitude) + for diagram in matrix_element.get('diagrams') + for amplitude in diagram.get('amplitudes')) + unrollable = diagram_generation.get_unrollable_quartic_vertices( + self.base_model) + + # every sum is a quartic current against the cubic pair it splits into + for cubic, quartic, coeff in sums: + self.assertTrue(helas_objects.HelasMatrixElement.is_unrolled_pair( + quartic, cubic, unrollable, self.cubic_ids)) + + # and every amplitude taken away is the target with some of the + # substitutions applied, which is exactly what reading the sums gives + for source in folded: + self.assertIn(source, merges) + target, coeff = merges[source] + self.assertIn(target, uses) + swap = dict((cubic, sums[index][1].get('number')) + for cubic, index in uses[target].items()) + mothers = [mother.get('number') + for mother in amplitudes[target].get('mothers')] + options = [] + for size in range(1, len(swap) + 1): + for combination in itertools.combinations(sorted(swap), size): + options.append(sorted(swap[number] + if number in combination else number + for number in mothers)) + self.assertIn(sorted(mother.get('number') for mother in + amplitudes[source].get('mothers')), options) + self.assertEqual(amplitudes[source].get('interaction_id'), + amplitudes[target].get('interaction_id')) + + # a folded amplitude must not also be a target + self.assertFalse(folded & set(uses)) + + def check_auto_order(self, initial, final): + """'auto' has to generate the 'speed' order, and 'slots' has to be + that same order reversed. + + This is what lets an output pick between the two: by then the diagrams + exist and only their order can still be changed, so the two modes are + only interchangeable at that point if one is the other backwards.""" + + def tags(mode): + madgraph.merge_quartic_vertices = mode + amplitude = diagram_generation.Amplitude(base_objects.Process( + {'legs':base_objects.LegList( + [base_objects.Leg({'id':pdg, 'state':False}) + for pdg in initial] + + [base_objects.Leg({'id':pdg, 'state':True}) + for pdg in final]), + 'model':self.base_model})) + return [str(diagram_generation.UnrollDiagramTag( + diagram, self.base_model, len(initial))) + for diagram in amplitude.get('diagrams')] + + speed, slots, auto = tags('speed'), tags('slots'), tags('auto') + self.assertEqual(slots, speed[::-1]) + # and it is a reordering, nothing gained or lost + self.assertEqual(sorted(slots), sorted(speed)) + self.assertEqual(len(set(speed)), len(speed)) + # 'auto' only takes that order once the process is big enough to pay + # for it, and otherwise leaves the generation order alone + if len(initial) + len(final) >= madgraph.merge_quartic_min_legs: + self.assertEqual(auto, speed) + else: + self.assertEqual(auto, tags(False)) + + def test_auto_order_gg_ggg(self): + self.check_auto_order([21, 21], [21, 21, 21]) + + def test_auto_order_gg_gggg(self): + self.check_auto_order([21, 21], [21, 21, 21, 21]) + + def test_auto_current_sums(self): + """'auto' keeps the sums, being the 'speed' order; 'slots' drops them + because reversing puts every target ahead of what feeds it""" + + import madgraph.core.helas_objects as helas_objects + + # five legs, so below merge_quartic_min_legs: 'auto' does not seed and + # so has no sums, while asking for 'speed' by name still does + for mode, wanted in (('auto', 0), ('speed', 7), ('slots', 0)): + madgraph.merge_quartic_vertices = mode + amplitude = diagram_generation.Amplitude(base_objects.Process( + {'legs':base_objects.LegList( + [base_objects.Leg({'id':21, 'state':False})] * 2 + + [base_objects.Leg({'id':21, 'state':True})] * 3), + 'model':self.base_model})) + element = helas_objects.HelasMatrixElement(amplitude) + self.assertEqual(len(element.get_quartic_current_sums()[0]), wanted) + + def test_auto_multiplicity_gate(self): + """'auto' leaves the seed rule off below merge_quartic_min_legs, where + it is measured not to pay for itself, and takes it above""" + + def seeded(nfinal, mode): + madgraph.merge_quartic_vertices = mode + amplitude = diagram_generation.Amplitude(base_objects.Process( + {'legs':base_objects.LegList( + [base_objects.Leg({'id':21, 'state':False})] * 2 + + [base_objects.Leg({'id':21, 'state':True})] * nfinal), + 'model':self.base_model})) + return bool(amplitude.seed_forbidden_cubic_ids) + + threshold = madgraph.merge_quartic_min_legs + for nfinal in (2, 3, 4): + wanted = (2 + nfinal) >= threshold + self.assertEqual(seeded(nfinal, 'auto'), wanted) + # asking for it by name is unconditional, that being the way to + # get the merging on a small process + self.assertTrue(seeded(nfinal, 'speed')) + self.assertTrue(seeded(nfinal, 'slots')) + self.assertFalse(seeded(nfinal, False)) + + def test_amplitude_slots(self): + """The AMP array is recycled: an amplitude summed into another frees + its entry, and the entries have to be reusable without two live + amplitudes ever sharing one""" + + import madgraph.core.helas_objects as helas_objects + + madgraph.merge_quartic_vertices = 'speed' + amplitude = diagram_generation.Amplitude(base_objects.Process( + {'legs':base_objects.LegList( + [base_objects.Leg({'id':21, 'state':False})] * 2 + + [base_objects.Leg({'id':21, 'state':True})] * 4), + 'model':self.base_model})) + element = helas_objects.HelasMatrixElement(amplitude) + slots, nslots, folds_at = element.get_amplitude_slots() + + self.assertEqual(element.get_number_of_amplitudes(), 510) + self.assertTrue(nslots < 510) + self.assertEqual(max(slots.values()), nslots) + self.assertEqual(min(slots.values()), 1) + + # replay the emission and check no entry is written while it still + # holds something with a reader to come + folded = set(element.get_quartic_current_sums()[2]) + order = [a.get('number') for d in element.get('diagrams') + for a in d.get('amplitudes') if a.get('number') not in folded] + self.assertEqual(len(order), len(slots)) + live = {} + for i, number in enumerate(order): + self.assertNotIn(slots[number], live) + live[slots[number]] = number + for target, source, coeff in folds_at.get(i, []): + self.assertIn(slots[source], live) + del live[slots[source]] + + def test_seed_inactive_by_default(self): + """Nothing changes unless madgraph.merge_quartic_vertices is set""" + + madgraph.merge_quartic_vertices = False + amplitude = diagram_generation.Amplitude(base_objects.Process( + {'legs':base_objects.LegList( + [base_objects.Leg({'id':21, 'state':False})] * 2 + + [base_objects.Leg({'id':21, 'state':True})] * 2), + 'model':self.base_model})) + self.assertEqual(amplitude.seed_forbidden_cubic_ids, frozenset()) + self.assertEqual(len(amplitude.get('diagrams')), 4) diff --git a/tests/unit_tests/fks/test_fks_base.py b/tests/unit_tests/fks/test_fks_base.py index f975840a3..9c6a39e92 100755 --- a/tests/unit_tests/fks/test_fks_base.py +++ b/tests/unit_tests/fks/test_fks_base.py @@ -23,6 +23,7 @@ import tests.unit_tests as unittest import madgraph.various.misc as misc +import madgraph import madgraph.fks.fks_base as fks_base import madgraph.fks.fks_common as fks_common import madgraph.core.base_objects as MG @@ -35,6 +36,15 @@ class TestFKSProcess(unittest.TestCase): """a class to test FKS Processes""" + def setUp(self): + # these build the fks amplitudes by hand, so they have to turn the + # four gluon merging off themselves, as FKSMultiProcess does + self.merge_quartic = madgraph.merge_quartic_vertices + madgraph.merge_quartic_vertices = False + + def tearDown(self): + madgraph.merge_quartic_vertices = self.merge_quartic + # the model, import the SM but remove 2nd and 3rd gen quarks remove_list = [3,4,5,6,-3,-4,-5,-6] mymodel = import_ufo.import_model('sm', options={'apply_flavor_grouping':False}) diff --git a/tests/unit_tests/fks/test_fks_common.py b/tests/unit_tests/fks/test_fks_common.py index 879f4d646..79919e196 100755 --- a/tests/unit_tests/fks/test_fks_common.py +++ b/tests/unit_tests/fks/test_fks_common.py @@ -24,6 +24,7 @@ sys.path.insert(0, os.path.join(root_path,'..','..')) import tests.unit_tests as unittest +import madgraph import madgraph.fks.fks_common as fks_common import madgraph.core.base_objects as MG import madgraph.core.color_algebra as color @@ -2933,9 +2934,16 @@ class TestLinkRBConfHEFT(unittest.TestCase): (only processes with 3 point interactions)""" def setUp(self): + # link_rb_configs reads the vertex decomposition of the real diagrams, + # so it runs with the four gluon merging off, as FKSMultiProcess does + self.merge_quartic = madgraph.merge_quartic_vertices + madgraph.merge_quartic_vertices = False if not hasattr(self, 'base_model'): TestLinkRBConfHEFT.base_model = import_ufo.import_model('heft') + def tearDown(self): + madgraph.merge_quartic_vertices = self.merge_quartic + def test_link_gghg_ggh(self): """tests that the real emission process gg>hg and born process gg>h are @@ -3062,9 +3070,16 @@ class TestLinkRBConfSM(unittest.TestCase): (only processes with 3 point interactions)""" def setUp(self): + # link_rb_configs reads the vertex decomposition of the real diagrams, + # so it runs with the four gluon merging off, as FKSMultiProcess does + self.merge_quartic = madgraph.merge_quartic_vertices + madgraph.merge_quartic_vertices = False if not hasattr(self, 'base_model'): TestLinkRBConfSM.base_model = import_ufo.import_model('sm', options={'apply_flavor_grouping':False}) + def tearDown(self): + madgraph.merge_quartic_vertices = self.merge_quartic + def test_link_udxwpg_udxwp(self): """tests that the real emission process ud~>w+g and born process u u~>w+ are correctly linked""" diff --git a/tests/unit_tests/loop/test_loop_helas_objects.py b/tests/unit_tests/loop/test_loop_helas_objects.py index 261b091e9..a2b6cf591 100755 --- a/tests/unit_tests/loop/test_loop_helas_objects.py +++ b/tests/unit_tests/loop/test_loop_helas_objects.py @@ -296,6 +296,11 @@ def check_HME_individual_diag_sanity(self,Amplitude, process,\ amp_number_apparition=[] for jamp in color_amplitudes: amp_number_apparition.extend([a[1] for a in jamp]) + # A four gluon contribution summed into another amplitude is + # dropped from the jamps on purpose, its colour factor being + # carried by the amplitude it was summed into + amp_number_apparition.extend( + myME.get_quartic_amplitude_merges().keys()) diagIndex=0 for i, diag in enumerate(diagSelection):