Implement Feistel-based permutation for LHS and convert samplers to stateless RNG

framework/include/samplers/Sampler.h

@@ -20,6 +20,7 @@
 #include "MultiApp.h"
 #include "VectorPostprocessorInterface.h"
 #include "ReporterInterface.h"
+#include "MooseRandomStateless.h"
 
 /**
  * This is the base class for Samplers as used within the Stochastic Tools module.
@@ -138,22 +139,6 @@ class Sampler : public MooseObject,
   libMesh::Parallel::Communicator & getLocalComm() { return _local_comm; }
 
 protected:
-  /**
-   * Enum describing the type of parallel communication to perform.
-   *
-   * Some routines require specific communication methods that not all processors
-   * see, these IDs will determine how that routine is performed:
-   *  - NONE routine is not distrubuted and things all can happen locally
-   *  - LOCAL routine is distributed on all processors
-   *  - SEMI_LOCAL routine is distributed only on processors that own rows
-   */
-  enum CommMethod
-  {
-    NONE = 0,
-    LOCAL = 1,
-    SEMI_LOCAL = 2
-  };
-
   // The following methods are the basic methods that should be utilized my most application
   // developers that are creating a custom Sampler.
 
@@ -175,24 +160,27 @@ class Sampler : public MooseObject,
   void setNumberOfRandomSeeds(std::size_t number);
 
   /**
-   * Get the next random number from the generator.
+   * Get nth random number from the generator.
+   * @param n 0-based index of the random number to generate
    * @param index The index of the seed, by default this is zero. To add additional seeds
    *              indices call the setNumberOfRequiedRandomSeeds method.
    *
    * @return A double for the random number, this is double because MooseRandom class uses double.
    */
-  Real getRand(unsigned int index = 0);
+  Real getRand(std::size_t n, unsigned int index = 0) const;
 
   /**
-   * Get the next random integer from the generator within the specified range [lower, upper)
-   * @param index The index of the seed, by default this is zero. To add additional seeds
-   *              indices call the setNumberOfRequiedRandomSeeds method.
+   * Get nth random integer from the generator within the specified range [lower, upper)
+   * @param n 0-based index of the random number to generate
    * @param lower Lower bounds
    * @param upper Upper bounds
+   * @param index The index of the seed, by default this is zero. To add additional seeds
+   *              indices call the setNumberOfRequiedRandomSeeds method.
    *
    * @return A integer for the random number
    */
-  uint32_t getRandl(unsigned int index, uint32_t lower, uint32_t upper);
+  unsigned int
+  getRandl(std::size_t n, unsigned int lower, unsigned int upper, unsigned int index = 0) const;
 
   /**
    * Base class must override this method to supply the sample distribution data.
@@ -202,17 +190,6 @@ class Sampler : public MooseObject,
    */
   virtual Real computeSample(dof_id_type row_index, dof_id_type col_index) = 0;
 
-  ///@{
-  /**
-   * Setup method called prior and after looping through distributions.
-   *
-   * These methods should not be called directly, each is automatically called by the public
-   * getGlobalSamples() or getLocalSamples() methods.
-   */
-  virtual void sampleSetUp(const SampleMode /*mode*/) {}
-  virtual void sampleTearDown(const SampleMode /*mode*/) {}
-  ///@}
-
   // The following methods are advanced methods that should not be needed by application developers,
   // but exist for special cases.
 
@@ -241,44 +218,23 @@ class Sampler : public MooseObject,
 
   /**
    * Method for advancing the random number generator(s) by the supplied number or calls to rand().
-   *
-   * TODO: This should be updated if the If the random number generator is updated to type that
-   * supports native advancing.
    */
   virtual void advanceGenerators(const dof_id_type count);
   virtual void advanceGenerator(const unsigned int seed_index, const dof_id_type count);
   void setAutoAdvanceGenerators(const bool state);
 
-  /**
-   * Helper for shuffling a vector of data in-place; the default assumes data is distributed
-   *
-   * NOTE: This will advance the generator by the size of the supplied vector.
-   */
-  template <typename T>
-  void shuffle(std::vector<T> & data,
-               const std::size_t seed_index = 0,
-               const CommMethod method = CommMethod::LOCAL);
-
   //@{
   /**
    * Callbacks for before and after execute.
    *
    * These were added to support of dynamic sampler sizes. Recall that execute is simply to advance
-   * the state of the generator such that the next sample will be unique. These methods allow
+   * the state of the generators such that the next sample will be unique. These methods allow
    * operations before and after the call to generator advancement.
    */
   virtual void executeSetUp() {}
   virtual void executeTearDown() {}
   ///@}
 
-  //@{
-  /**
-   * Here we save/restore generator states
-   */
-  void saveGeneratorState() { _generator.saveState(); }
-  void restoreGeneratorState() { _generator.restoreState(); }
-  //@}
-
   /**
    * This is where the sampler partitioning is defined. It is NOT recommended to
    * override this function unless you know EXACTLY what you are doing
@@ -311,11 +267,8 @@ class Sampler : public MooseObject,
                                         const std::string & name,
                                         InputParameters & parameters);
   /**
-   * Store the state of the MooseRandom generator so that new calls to
-   * getGlobalSamples/getLocalSamples methods will create new numbers.
-   *
-   * The execute() method is called in the init() method of this class and
-   * FEProblemBase::executeSamplers; it should not be called elsewhere.
+   * Advance MooseRandomStateless generators so that new calls to
+   * sample methods will create new numbers.
    */
   void execute();
   friend void FEProblemBase::objectExecuteHelper<Sampler>(const std::vector<Sampler *> & objects);
@@ -330,8 +283,8 @@ class Sampler : public MooseObject,
    */
   void advanceGeneratorsInternal(const dof_id_type count);
 
-  /// Random number generator, don't give users access. Control it via the interface from this class.
-  MooseRandom _generator;
+  /// Random number generators, don't give users access. Control it via the interface from this class.
+  std::vector<std::unique_ptr<MooseRandomStateless>> _generators;
 
   /// Number of rows for this processor
   dof_id_type _n_local_rows;
@@ -382,15 +335,3 @@ class Sampler : public MooseObject,
   /// Flag for disabling automatic generator advancing
   bool _auto_advance_generators;
 };
-
-template <typename T>
-void
-Sampler::shuffle(std::vector<T> & data, const std::size_t seed_index, const CommMethod method)
-{
-  if (method == CommMethod::NONE)
-    MooseUtils::shuffle<T>(data, _generator, seed_index, nullptr);
-  else if (method == CommMethod::LOCAL)
-    MooseUtils::shuffle<T>(data, _generator, seed_index, &_communicator);
-  else if (method == CommMethod::SEMI_LOCAL)
-    MooseUtils::shuffle<T>(data, _generator, seed_index, &_local_comm);
-}

framework/include/utils/MooseRandomPerturbation.h

@@ -0,0 +1,127 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "MooseError.h"
+
+#include <cstdint>
+
+/**
+ * Generates a keyed pseudo-random permutation of the integers [0, n) using a
+ * balanced Feistel network. Given the same seed and n, the mapping is fully
+ * deterministic and bijective: every input in [0, n) maps to a unique output
+ * in [0, n). Different seeds produce statistically independent permutations.
+ *
+ * Because n need not be a power of two, the Feistel network operates on a
+ * padded domain of size 2^(2*half_bits) >= n and uses cycle-walking: if the
+ * raw output falls outside [0, n) it is re-applied until a valid index is
+ * reached. This guarantees termination because the padded permutation is
+ * itself a bijection and n > 0 ensures at least one valid output exists.
+ *
+ * The permutation is also invertible via invert(), which runs the Feistel
+ * rounds in reverse order.
+ */
+class MooseRandomPerturbation
+{
+public:
+  /**
+   * Construct a permutation of [0, \p n) keyed by \p seed.
+   *
+   * @param seed  64-bit seed; split into two 32-bit subkeys for the round function.
+   * @param n     Size of the domain to permute; must be > 0.
+   * @param rounds Number of Feistel rounds. More rounds improve mixing at the
+   *               cost of throughput; 8 is sufficient for sampling applications.
+   */
+  MooseRandomPerturbation(uint64_t seed, unsigned int n, unsigned int rounds = 8);
+
+  /**
+   * Map \p x to its permuted index in [0, n).
+   *
+   * Applies the Feistel network to the padded domain and cycle-walks until
+   * the result falls within [0, n).
+   *
+   * @param x  Input index; must satisfy x < n.
+   * @return   A unique index in [0, n). Calling permute for every x in [0, n)
+   *           yields each value in [0, n) exactly once.
+   */
+  uint32_t permute(uint32_t x) const;
+
+  /**
+   * Recover the original index from a permuted value, i.e. invert(permute(x)) == x.
+   *
+   * Runs the Feistel rounds in reverse order and cycle-walks back into [0, n).
+   *
+   * @param y  Permuted index; must satisfy y < n.
+   * @return   The unique x in [0, n) such that permute(x) == y.
+   */
+  uint32_t invert(uint32_t y) const;
+
+private:
+  /**
+   * Apply one full pass of the balanced Feistel network over the padded
+   * domain [0, 2^(2*half_bits)). The result may fall outside [0, n); the
+   * caller is responsible for cycle-walking.
+   *
+   * @param x  Input value in the padded domain.
+   * @return   Permuted value in the padded domain.
+   */
+  uint32_t permutePadded(uint32_t x) const;
+
+  /**
+   * Invert one full pass of the Feistel network by running rounds in reverse.
+   *
+   * @param y  Value in the padded domain produced by permutePadded.
+   * @return   The original input to permutePadded that produced \p y.
+   */
+  uint32_t invertPadded(uint32_t y) const;
+
+  /**
+   * Keyed round function F(half, round) used in each Feistel step.
+   *
+   * Mixes the half-block with both subkeys and a round-dependent constant
+   * derived from the golden-ratio fractional bits (0x9e3779b9), then passes
+   * the result through a 32-bit avalanche hash (mix32).
+   *
+   * @param half   The half-block value (R in the forward direction).
+   * @param round  Zero-based round index, used to vary the constant each round.
+   * @return       Mixed value masked to \p _half_bits width.
+   */
+  uint32_t roundFunction(uint32_t half, unsigned int round) const;
+
+  /**
+   * Return the number of bits needed to represent values in [0, n-1], i.e.
+   * ceil(log2(n)). Returns 0 for n == 1 (no bits needed to index a single element).
+   *
+   * @param n  Must be > 0.
+   */
+  static unsigned int ceilLog2(uint32_t n);
+
+  /**
+   * Bijective 32-bit avalanche hash (finalizer from Murmur3 / degski hash).
+   * Used inside the round function to achieve good bit diffusion.
+   *
+   * @param x  32-bit input.
+   * @return   Hashed 32-bit output; the mapping is invertible.
+   */
+  static uint32_t mix32(uint32_t x);
+
+  /// Lower 32 bits of the seed, used as the first subkey in the round function
+  const uint32_t _k0;
+  /// Upper 32 bits of the seed, used as the second subkey in the round function
+  const uint32_t _k1;
+  /// Size of the permutation domain [0, n)
+  const unsigned int _n;
+  /// Number of bits in each Feistel half-block: ceil((ceil(log2(n)) + 1) / 2)
+  const unsigned int _half_bits;
+  /// Bitmask of width _half_bits, used to keep half-block arithmetic in range
+  const uint32_t _half_mask;
+  /// Number of Feistel rounds to apply per permute/invert call
+  const unsigned int _rounds;
+};

framework/include/utils/MooseRandomStateless.h

@@ -146,7 +146,7 @@ class MooseRandomStateless
     const auto range = upper - lower;
 
     auto [it, is_new] = _randlb_generators.try_emplace(
-        range, [&range](mt_state * state) { return rds_iuniform(state, 0, range); }, _seed);
+        range, [range](mt_state * state) { return rds_iuniform(state, 0, range); }, _seed);
 
     if (is_new)
       it->second.advance(_advance_count);

framework/src/problems/FEProblemBase.C

@@ -5035,13 +5035,7 @@ FEProblemBase::execute(const ExecFlagType & exec_type)
         exec_type == EXEC_SUBDOMAIN || exec_type == EXEC_NONLINEAR || exec_type == EXEC_LINEAR))
     customSetup(exec_type);
 
-  // Samplers; EXEC_INITIAL is not called because the Sampler::init() method that is called after
-  // construction makes the first Sampler::execute() call. This ensures that the random number
-  // generator object is the correct state prior to any other object (e.g., Transfers) attempts to
-  // extract data from the Sampler. That is, if the Sampler::execute() call is delayed to here
-  // then it is not in the correct state for other objects.
-  if (exec_type != EXEC_INITIAL)
-    executeSamplers(exec_type);
+  executeSamplers(exec_type);
 
   // Pre-aux UserObjects
   computeUserObjects(exec_type, Moose::PRE_AUX);