Update: optional antithetic pairs of random numbers

Author: heptaflar

README.md

@@ -24,6 +24,7 @@ will print a detailed usage messages to the screen, like
     usage: kaleidoscope-scatter [-h] --source-type
                                 {esa-cci-oc,esa-scope-cs,esa-scope-pp,ghrsst,glorys}
                                 --selector {0,1,2,...,100}
+                                [--antithetic]
                                 [--engine-reader {h5netcdf,netcdf4,zarr}]
                                 [--engine-writer {h5netcdf,netcdf4,zarr}]
                                 [--log-level {debug,info,warning,error,off}]
@@ -46,6 +47,8 @@ will print a detailed usage messages to the screen, like
       --selector {0,1,2,...,100}
                             the Monte Carlo stream selector. An integral number
                             which must not be negative. (default: None)
+      --antithetic          enable pairwise antithetic Monte Carlo simulation.
+                            (default: False)
       --engine-reader {h5netcdf,netcdf4,zarr}
                             specify the engine used to read the source product
                             file. (default: None)

kaleidoscope/config/config.yml

@@ -1,8 +1,8 @@
 #  Copyright (c) Brockmann Consult GmbH, 2025
 #  License: MIT
 
-## The default source type.
-source_type: _
+## No pairwise antithetic Monte Carlo simulation.
+antithetic: False
 
 ## Detect the reader engine automatically
 engine_reader:

kaleidoscope/main/scatter.py

@@ -111,6 +111,13 @@ def _add_options(parser):
             required=True,
             dest="selector",
         )
+        parser.add_argument(
+            "--antithetic",
+            help="enable pairwise antithetic Monte Carlo simulation.",
+            action="store_true",
+            required=False,
+            dest="antithetic",
+        )
         parser.add_argument(
             "--engine-reader",
             help="specify the engine used to read the source product file.",

kaleidoscope/operators/scatterop.py

@@ -108,7 +108,12 @@ def randomize(
                 z = da.clip(z, config["clip"][0], config["clip"][1])
         elif "uncertainty" in config:
             s = self.seed(self.uuid(v))
-            f = Randomize(m=x.ndim, dist=config["distribution"], seed=s)
+            f = Randomize(
+                m=x.ndim,
+                dist=config["distribution"],
+                seed=s,
+                antithetic=self._args.antithetic,
+            )
             u = (
                 target[config["uncertainty"]]
                 if isinstance(config["uncertainty"], str)
@@ -130,7 +135,12 @@ def randomize(
             )
         else:
             s = self.seed(self.uuid(v))
-            f = Randomize(m=x.ndim, dist=config["distribution"], seed=s)
+            f = Randomize(
+                m=x.ndim,
+                dist=config["distribution"],
+                seed=s,
+                antithetic=self._args.antithetic,
+            )
             b = target[config["bias"]]
             r = target[config["rmsd"]]
             z = f.apply_to(
@@ -177,7 +187,11 @@ def seed(self, uuid: uuid.UUID, n: int = 4) -> np.ndarray:
         """
         from numpy.random import Philox
 
-        seed = uuid.int + self._args.selector
+        seed = uuid.int + (
+            self._args.selector
+            if not self._args.antithetic
+            else (self._args.selector + 1) // 2
+        )
         g = DefaultGenerator(Philox(seed))
         return np.array([g.next() for _ in range(n)], dtype=np.int64)