Implement reproducible event generation

.github/workflows/acceptancetest_mg7.yml

@@ -201,6 +201,92 @@ jobs:
             export PATH="$HOME/.cache/HEPtools/bin:$PATH"
             ./tests/test_manager.py test_group_subprocess_mg7 -pA -t0 -l INFO
 
+  acceptancetest_mg7_vegas_reproducibility:
+    needs: build_madspace
+    # mg7 seeded-generation reproducibility (p p > t t~, plain VEGAS-optimized
+    # run): the same run_card seed must give a byte-identical LHE file across
+    # independent runs, a different seed must not. Self-skips if the madspace
+    # + LHAPDF(NNPDF23) stack is unavailable.
+    runs-on: ubuntu-24.04
+    if: github.event_name == 'push' || github.event.pull_request.head.repo.fork == true
+    steps:
+      - uses: actions/checkout@v4
+      - uses: ./.github/actions/checkout_mg5
+      - uses: ./.github/actions/install_madspace
+      - uses: ./.github/actions/restore-pip-cache
+      - uses: ./.github/actions/restore_heptools
+      - name: test one of the test test_vegas_reproducibility_mg7
+        run: |
+            cd $GITHUB_WORKSPACE
+            export PATH="$HOME/.cache/HEPtools/bin:$PATH"
+            ./tests/test_manager.py test_vegas_reproducibility_mg7 -pA -t0 -l INFO
+
+  acceptancetest_mg7_madnis_reproducible_mode:
+    needs: build_madspace
+    # mg7 madnis.reproducible = True (p p > t t~): madnis training itself must be
+    # byte-identical for the same seed independent of the cpu thread pool size,
+    # both with online-only training and with buffered (off-policy replay)
+    # training enabled. Self-skips if the madspace + LHAPDF(NNPDF23) stack is
+    # unavailable.
+    runs-on: ubuntu-24.04
+    if: github.event_name == 'push' || github.event.pull_request.head.repo.fork == true
+    steps:
+      - uses: actions/checkout@v4
+      - uses: ./.github/actions/checkout_mg5
+      - uses: ./.github/actions/install_madspace
+      - uses: ./.github/actions/restore-pip-cache
+      - uses: ./.github/actions/restore_heptools
+      - name: test one of the test test_madnis_reproducible_mode_mg7
+        run: |
+            cd $GITHUB_WORKSPACE
+            export PATH="$HOME/.cache/HEPtools/bin:$PATH"
+            ./tests/test_manager.py test_madnis_reproducible_mode_mg7 -pA -t0 -l INFO
+
+  acceptancetest_mg7_gridpack_reproducibility:
+    needs: build_madspace
+    # mg7 seeded-generation reproducibility for a gridpack trained with
+    # madnis (p p > t t~): the gridpack's own bin/generate_events --seed must
+    # give a byte-identical LHE file across independent runs with the same
+    # seed, and a different result for a different seed. Training itself is
+    # not required to be deterministic. Self-skips if the madspace +
+    # LHAPDF(NNPDF23) stack is unavailable.
+    runs-on: ubuntu-24.04
+    if: github.event_name == 'push' || github.event.pull_request.head.repo.fork == true
+    steps:
+      - uses: actions/checkout@v4
+      - uses: ./.github/actions/checkout_mg5
+      - uses: ./.github/actions/install_madspace
+      - uses: ./.github/actions/restore-pip-cache
+      - uses: ./.github/actions/restore_heptools
+      - name: test one of the test test_gridpack_reproducibility_mg7
+        run: |
+            cd $GITHUB_WORKSPACE
+            export PATH="$HOME/.cache/HEPtools/bin:$PATH"
+            ./tests/test_manager.py test_gridpack_reproducibility_mg7 -pA -t0 -l INFO
+
+  acceptancetest_mg7_gridpack_reproducibility_vegas:
+    needs: build_madspace
+    # mg7 seeded-generation reproducibility for a plain VEGAS-optimized
+    # gridpack (p p > t t~, no madnis training): the gridpack's own
+    # bin/generate_events --seed must give a byte-identical LHE file across
+    # independent runs with the same seed, and a different result for a
+    # different seed. Companion to acceptancetest_mg7_gridpack_reproducibility,
+    # which covers the madnis-trained case. Self-skips if the madspace +
+    # LHAPDF(NNPDF23) stack is unavailable.
+    runs-on: ubuntu-24.04
+    if: github.event_name == 'push' || github.event.pull_request.head.repo.fork == true
+    steps:
+      - uses: actions/checkout@v4
+      - uses: ./.github/actions/checkout_mg5
+      - uses: ./.github/actions/install_madspace
+      - uses: ./.github/actions/restore-pip-cache
+      - uses: ./.github/actions/restore_heptools
+      - name: test one of the test test_gridpack_reproducibility_vegas_mg7
+        run: |
+            cd $GITHUB_WORKSPACE
+            export PATH="$HOME/.cache/HEPtools/bin:$PATH"
+            ./tests/test_manager.py test_gridpack_reproducibility_vegas_mg7 -pA -t0 -l INFO
+
   acceptancetest_mg7_merged_flavor_uq:
     needs: build_madspace
     # mg7 cross-section for the merged-flavor u q > u q (q = u d), pinned to the

madgraph/iolibs/template_files/mg7/gridpack.py

@@ -31,7 +31,6 @@
 import tomllib
 import argparse
 
-
 def resolve_verbosity(verbosity: str) -> str:
     """Resolve the run_card "auto" verbosity to "pretty"/"log" depending on
     whether stdout is attached to a terminal; other values pass through
@@ -69,6 +68,12 @@ def main() -> None:
     # parse command line arguments
     parser = argparse.ArgumentParser()
     parser.add_argument("--run_name", type=str, default=run_args["run_name"])
+    parser.add_argument(
+        "--seed", type=int, default=run_args.get("seed", 0),
+        help="0 draws a fresh random seed; a positive integer makes the run "
+             "reproducible at any --cpu_thread_pool_size (for --output_format lhe "
+             "also set combine_thread_pool_size = 1)"
+    )
     parser.add_argument("--device", type=str, nargs="*")
     parser.add_argument(
         "--cpu_thread_pool_size", type=int, default=run_args["cpu_thread_pool_size"]
@@ -173,6 +178,7 @@ def main() -> None:
         channels=channel_generators,
         status_file=os.path.join(run_path, "info.json"),
         config=config,
+        seed=args.seed,
     )
 
     # run generation

madgraph/iolibs/template_files/mg7/madevent.py

@@ -212,7 +212,7 @@ def init_context(self) -> None:
         self.device_types = []
         self.devices = []
         self.pool_sizes = []
-        for i, device_name in enumerate(device_names):
+        for device_name in device_names:
             if ":" in device_name:
                 device_type, device_index_str = device_name.split(":")
                 device_index = int(device_index_str)
@@ -426,6 +426,9 @@ def build_event_generator(self, phasespaces: list[PhaseSpace]) -> ms.EventGenera
             channels=channel_generators,
             status_file=os.path.join(self.run_path, "info.json"),
             config=self.event_generator_config,
+            # 0 conventionally means "no seed" in the run card; EventGenerator takes
+            # an actual seed or None (non-deterministic), not a 0 sentinel.
+            seed=self.run_card["run"]["seed"] or None,
         )
         unused_globals = (
             set(self.contexts[0].global_names()) - event_generator.used_globals()
@@ -436,44 +439,50 @@ def build_event_generator(self, phasespaces: list[PhaseSpace]) -> ms.EventGenera
         return event_generator
 
     def survey_phasespaces(
-        self, phasespaces: list[PhaseSpace | None]
+        self, phasespaces: list[PhaseSpace | None], survey_pass: int = 0
     ) -> ms.EventGenerator | None:
         ps_filtered = [ps for ps in phasespaces if ps is not None]
         if len(ps_filtered) == 0:
             return None
         event_generator = self.build_event_generator(ps_filtered)
-        event_generator.survey()
+        event_generator.survey(survey_pass)
         return event_generator
 
     def survey(self) -> None:
+        # survey_pass distinguishes the survey() calls below: "both" mode can
+        # re-survey a channel carried over unchanged from the multichannel pass
+        # into the final (simplified) pass, and both passes schedule jobs on the
+        # same underlying ChannelEventGenerator. The explicit pass index keeps each
+        # pass's job seeds independent of the other passes' job counts, rather than
+        # depending on call history.
         phasespace_mode = self.run_card["phasespace"]["mode"]
         if phasespace_mode == "multichannel":
             self.phasespaces = [
                 subproc.build_multichannel_phasespace()
                 for subproc in self.subprocesses
             ]
-            self.event_generator = self.survey_phasespaces(self.phasespaces)
+            self.event_generator = self.survey_phasespaces(self.phasespaces, 0)
         elif phasespace_mode == "flat":
             self.phasespaces = [
                 subproc.build_flat_phasespace()
                 for subproc in self.subprocesses
             ]
-            self.event_generator = self.survey_phasespaces(self.phasespaces)
+            self.event_generator = self.survey_phasespaces(self.phasespaces, 0)
         elif phasespace_mode == "both":
             kept_count = self.run_card["phasespace"]["simplified_channel_count"]
             phasespaces_multi = [
                 subproc.build_multichannel_phasespace()
                 for subproc in self.subprocesses
             ]
-            evgen_multi = self.survey_phasespaces(phasespaces_multi)
+            evgen_multi = self.survey_phasespaces(phasespaces_multi, 0)
 
             phasespaces_flat = [
                 subproc.build_flat_phasespace()
                 if len(subproc.meta["channels"]) > kept_count + 1 else
                 None
                 for subproc in self.subprocesses
             ]
-            #evgen_flat = self.survey_phasespaces(phasespaces_flat, "flat")
+            #evgen_flat = self.survey_phasespaces(phasespaces_flat, 1)
 
             channel_status = evgen_multi.channel_status()
             cross_sections = []
@@ -494,7 +503,7 @@ def survey(self) -> None:
                     self.subprocesses, phasespaces_multi, phasespaces_flat, cross_sections
                 )
             ]
-            self.event_generator = self.survey_phasespaces(self.phasespaces)
+            self.event_generator = self.survey_phasespaces(self.phasespaces, 2)
         else:
             raise ValueError("Unknown phasespace mode")
 
@@ -534,6 +543,7 @@ def train_madnis(self) -> None:
         config.buffer_unweighting_quantile = madnis_args["buffer_unweighting_quantile"]
         config.fixed_cwnet_fraction = madnis_args["fixed_cwnet_fraction"]
         config.softclip_threshold = madnis_args["softclip_threshold"]
+        config.reproducible = madnis_args["reproducible"]
         madnis_phasespaces = []
         integrands = []
         cwnets = []
@@ -559,6 +569,11 @@ def train_madnis(self) -> None:
             config=config,
             integrands=integrands,
             cwnets=cwnets,
+            # Reuses the run_card's own seed (also used by build_event_generator()).
+            # Only the single-channel CPU sample-generation path is currently seeded
+            # -- buffered training and GPU multi-channel batches are still
+            # non-deterministic.
+            seed=run_args["seed"] or None,
         )
         madnis_training.train()
         for phasespace, active_channels in zip(
@@ -1301,6 +1316,9 @@ def simplify_phasespace(
 
     def build_madnis(self, phasespace: PhaseSpace) -> PhaseSpace:
         madnis_args = self.process.run_card["madnis"]
+        # Shared across all networks below: initialize_globals() derives an
+        # independent, non-colliding stream per tensor from this one base seed.
+        seed = self.process.run_card["run"]["seed"] or None
         channels = []
         for channel_id, channel in enumerate(phasespace.channels):
             prefix = f"subproc{self.subproc_id}.channel{channel_id}"
@@ -1318,7 +1336,7 @@ def build_madnis(self, phasespace: PhaseSpace) -> PhaseSpace:
                     subnet_layers=madnis_args["discrete_layers"],
                     subnet_activation=self.activation(madnis_args["discrete_activation"]),
                 )
-                discrete_before.initialize_globals(self.process.contexts[0])
+                discrete_before.initialize_globals(self.process.contexts[0], seed)
                 cond_dim += perm_count
 
             flow_dim = channel.phasespace_mapping.random_dim()
@@ -1333,10 +1351,11 @@ def build_madnis(self, phasespace: PhaseSpace) -> PhaseSpace:
                 invert_spline=madnis_args["flow_invert_spline"],
             )
             if channel.adaptive_mapping is None:
-                flow.initialize_globals(self.process.contexts[0])
+                flow.initialize_globals(self.process.contexts[0], seed)
             else:
                 flow.initialize_from_vegas(
-                    self.process.contexts[0], channel.adaptive_mapping.grid_name()
+                    self.process.contexts[0], channel.adaptive_mapping.grid_name(),
+                    seed
                 )
             cond_dim += flow_dim
 
@@ -1351,7 +1370,7 @@ def build_madnis(self, phasespace: PhaseSpace) -> PhaseSpace:
                     subnet_layers=madnis_args["discrete_layers"],
                     subnet_activation=self.activation(madnis_args["discrete_activation"]),
                 )
-                discrete_after.initialize_globals(self.process.contexts[0])
+                discrete_after.initialize_globals(self.process.contexts[0], seed)
 
             channels.append(Channel(
                 phasespace_mapping = channel.phasespace_mapping,
@@ -1420,7 +1439,9 @@ def build_cwnet(self, channel_count: int) -> ms.ChannelWeightNetwork:
             activation=self.activation(madnis_args["cwnet_activation"]),
             prefix=f"subproc{self.subproc_id}.cwnet",
         )
-        cwnet.initialize_globals(self.process.contexts[0])
+        cwnet.initialize_globals(
+            self.process.contexts[0], self.process.run_card["run"]["seed"] or None
+        )
         return cwnet
 
     def t_channel_mode(self, name: str) -> ms.PhaseSpaceMapping.TChannelMode:

madgraph/iolibs/template_files/mg7/run_card.toml

@@ -1,5 +1,11 @@
 [run]
 run_name = %(run.run_name)s
+# 0 draws a fresh random seed each run; a positive integer makes the run
+# reproducible: the same (seed, thread-pool sizes) reproduces bit-identically and
+# distinct seeds are statistically independent. Generation is reproducible at any
+# cpu_thread_pool_size; for output_format = lhe also set combine_thread_pool_size = 1
+# (compact_npy and lhe_npy are reproducible at any combine size).
+seed = %(run.seed)s
 # options: cuda, hip, cpp, cppnone, cppsse4, cppavx2, cpp512y, cpp512z, cppauto
 devices = %(run.devices)s
 # options:
@@ -136,3 +142,4 @@ batch_size_threshold = %(madnis.batch_size_threshold)s
 channel_grouping_mode = %(madnis.channel_grouping_mode)s # options: none, uniform, learned
 fixed_cwnet_fraction = %(madnis.fixed_cwnet_fraction)s
 softclip_threshold = %(madnis.softclip_threshold)s
+reproducible = %(madnis.reproducible)s # make training fully reproducible for a given seed, independent of thread count (CPU only, some throughput/memory cost)

madgraph/various/banner.py

@@ -6453,6 +6453,12 @@ def default_setup(self):
 
         # ----------------------------- [run] --------------------------
         self.add_toml_param('run', 'run_name', "run", gridpack=True)
+        self.add_toml_param('run', 'seed', 0, gridpack=True,
+            comment="0 draws a fresh random seed each run; any positive integer "
+                    "makes the run reproducible: the same (seed, thread-pool sizes) "
+                    "reproduces bit-identically. Generation is reproducible at any "
+                    "cpu_thread_pool_size; for output_format = lhe also set "
+                    "combine_thread_pool_size = 1")
         self.add_toml_param('run', 'devices', ["cppnone"], typelist=str, gridpack=True,
             comment="options: cuda, hip, cpp, cppnone, cppsse4, cppavx2, cpp512y, cpp512z, cppauto")
         self.add_toml_param('run', 'simd_vector_size', -1,
@@ -6588,6 +6594,7 @@ def default_setup(self):
             allowed=['none', 'uniform', 'learned'])
         self.add_toml_param('madnis', 'fixed_cwnet_fraction', 0.33)
         self.add_toml_param('madnis', 'softclip_threshold', 30.0)
+        self.add_toml_param('madnis', 'reproducible', False)
 
         # ----------------- dynamic (free-form) sections ---------------
         self.dynamic_sections['multiparticles'] = collections.OrderedDict([

madspace/.clang-format-ignore

@@ -7,3 +7,5 @@ src/cpu/runtime_mixin.inc
 src/gpu/runtime_mixin.inc
 src/cpu/runtime_backward_mixin.inc
 src/gpu/runtime_backward_mixin.inc
+include/madspace/mixmax/mixmax.hpp
+include/madspace/mixmax/mixmax_skip_N17.c

madspace/CMakeLists.txt

@@ -180,6 +180,7 @@ add_library(
     src/driver/backend.cpp
     src/driver/lhe_output.cpp
     src/driver/madnis_training.cpp
+    src/driver/random.cpp
     src/compgraphs/type.cpp
     src/compgraphs/function.cpp
     src/compgraphs/instruction.cpp
@@ -232,6 +233,7 @@ add_library(
     include/madspace/driver/lhe_output.hpp
     include/madspace/driver/logger.hpp
     include/madspace/driver/madnis_training.hpp
+    include/madspace/driver/random.hpp
     include/madspace/compgraphs.hpp
     include/madspace/compgraphs/type.hpp
     include/madspace/compgraphs/function.hpp

madspace/include/madspace/compgraphs/instruction.hpp

@@ -28,6 +28,9 @@ class Instruction {
     const std::string& name() const { return _name; }
     int opcode() const { return _opcode; }
     bool differentiable() const { return _differentiable; }
+    // draws from a shared RNG stream: never CSE'd, always scheduled on the main
+    // GPU stream (see FunctionBuilder::instruction)
+    virtual bool is_random() const { return false; }
 
 protected:
     void check_arg_count(const ValueVec& args, std::size_t count) const;
@@ -197,20 +200,23 @@ class RandomInstruction : public Instruction {
     RandomInstruction(int opcode, bool differentiable) :
         Instruction("random", opcode, differentiable) {}
     TypeVec signature(const ValueVec& args) const override;
+    bool is_random() const override { return true; }
 };
 
 class RandomIntInstruction : public Instruction {
 public:
     RandomIntInstruction(int opcode, bool differentiable) :
         Instruction("random_int", opcode, differentiable) {}
     TypeVec signature(const ValueVec& args) const override;
+    bool is_random() const override { return true; }
 };
 
 class UnweightInstruction : public Instruction {
 public:
     UnweightInstruction(int opcode, bool differentiable) :
         Instruction("unweight", opcode, differentiable) {}
     TypeVec signature(const ValueVec& args) const override;
+    bool is_random() const override { return true; }
 };
 
 class MatrixElementInstruction : public Instruction {