USACE-RMC · HadenSmith · Apr 6, 2026 · Apr 6, 2026
diff --git a/.github/workflows/NuGetPublish.yml b/.github/workflows/NuGetPublish.yml
@@ -7,6 +7,7 @@ on:
 jobs:
   publish:
     runs-on: windows-latest
+    timeout-minutes: 60
     steps:
       - uses: actions/checkout@v4
         with:
@@ -27,9 +28,6 @@ jobs:
       - name: Build
         run: dotnet build Numerics/Numerics.csproj -c Release -p:Version=${{ steps.version.outputs.VERSION }}
 
-      - name: Test (all frameworks)
-        run: dotnet test -c Release --no-build
-
       - name: Pack
         run: dotnet pack Numerics/Numerics.csproj -c Release --no-build -p:Version=${{ steps.version.outputs.VERSION }}
 

diff --git a/.github/workflows/Release.yml b/.github/workflows/Release.yml
@@ -9,7 +9,7 @@ jobs:
   build:
     uses: HydrologicEngineeringCenter/dotnet-workflows/.github/workflows/release.yml@main
     with:
-      dotnet-version: '9.0.x'
+      dotnet-version: '10.0.x'
       project-names: 'Numerics'
       run-tests: true
       nuget-source: 'https://www.hec.usace.army.mil/nexus/repository/consequences-nuget-public/'

diff --git a/.github/workflows/Snapshot.yml b/.github/workflows/Snapshot.yml
@@ -8,9 +8,9 @@ jobs:
   build:
     uses: HydrologicEngineeringCenter/dotnet-workflows/.github/workflows/snapshot.yml@main
     with:
-      dotnet-version: '9.0.x'
+      dotnet-version: '10.0.x'
       project-names: 'Numerics'
-      version: '1.0.0'
+      version: '2.0.0'
       run-tests: true
       nuget-source: 'https://www.hec.usace.army.mil/nexus/repository/consequences-nuget-public/'
     secrets:

diff --git a/Numerics/Distributions/Univariate/VonMises.cs b/Numerics/Distributions/Univariate/VonMises.cs
@@ -34,6 +34,7 @@
 using Numerics.Mathematics.Optimization;
 using Numerics.Mathematics.RootFinding;
 using Numerics.Mathematics.SpecialFunctions;
+using Numerics.Sampling;
 
 namespace Numerics.Distributions
 {
@@ -439,7 +440,7 @@ public override double[] GenerateRandomValues(int sampleSize, int seed = -1)
             if (_parametersValid == false)
                 ValidateParameters([Mu, Kappa], true);
 
-            var rng = seed < 0 ? new Random() : new Random(seed);
+            var rng = seed < 0 ? new MersenneTwister() : new MersenneTwister(seed);
             var values = new double[sampleSize];
 
             if (_kappa < 1e-10)

diff --git a/Numerics/Machine Learning/Supervised/KNearestNeighbors.cs b/Numerics/Machine Learning/Supervised/KNearestNeighbors.cs
@@ -30,6 +30,7 @@
 
 using Numerics.Data.Statistics;
 using Numerics.Mathematics.LinearAlgebra;
+using Numerics.Sampling;
 using System;
 using System.Linq;
 using System.Threading.Tasks;
@@ -363,7 +364,7 @@ public KNearestNeighbors(Matrix x, Vector y, int k)
         /// <param name="seed">Optional. The prng seed. If negative or zero, then the computer clock is used as a seed.</param>
         private double[]? kNNBootstrapPredict(Matrix xTrain, Vector yTrain, Matrix xTest, int seed = -1)
         {
-            var rnd = seed > 0 ? new Random(seed) : new Random();
+            var rnd = seed > 0 ? new MersenneTwister(seed) : new MersenneTwister();
             var idxs = rnd.NextIntegers(0, xTrain.NumberOfRows, xTrain.NumberOfRows);
             var bootX = new Matrix(xTrain.NumberOfRows, xTrain.NumberOfColumns);
             var bootY = new Vector(yTrain.Length);
@@ -393,7 +394,7 @@ public KNearestNeighbors(Matrix x, Vector y, int k)
             var output = new double[xTest.NumberOfRows, 4]; // lower, median, upper, mean
 
             var bootResults = new double[xTest.NumberOfRows, realizations];
-            var rnd = seed > 0 ? new Random(seed) : new Random();
+            var rnd = seed > 0 ? new MersenneTwister(seed) : new MersenneTwister();
             var seeds = rnd.NextIntegers(realizations);
 
             // Bootstrap the predictions

diff --git a/Numerics/Mathematics/Integration/MonteCarloIntegration.cs b/Numerics/Mathematics/Integration/MonteCarloIntegration.cs
@@ -28,6 +28,7 @@
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
 
+using Numerics.Sampling;
 using System;
 using System.Collections.Generic;
 using System.Linq;
@@ -83,7 +84,7 @@ public MonteCarloIntegration(Func<double[], double> function, int dimensions, IL
             Dimensions = dimensions;
             Min = min.ToArray();
             Max = max.ToArray();
-            Random = new Random();
+            Random = new MersenneTwister();
         }
 
         /// <summary>

diff --git a/Numerics/Mathematics/Optimization/Global/ShuffledComplexEvolution.cs b/Numerics/Mathematics/Optimization/Global/ShuffledComplexEvolution.cs
@@ -175,9 +175,9 @@ protected override void Optimize()
 
             // Initialize shuffled complex evolution (SCE) algorithm
             // Create list of complex prngs
-            var complexPRNGs = new List<Random>();
+            var complexPRNGs = new List<MersenneTwister>();
             for (i = 0; i < Complexes; i++)
-                complexPRNGs.Add(new Random(prng.Next()));
+                complexPRNGs.Add(new MersenneTwister(prng.Next()));
 
             // Create trapezoidal cumulative probability for points in complex
             var cdf = Trapezoidal(CCEIterations);
@@ -211,7 +211,7 @@ protected override void Optimize()
                 }
 
                 // Rank the points in terms of fitness (i.e., order of increasing function value).
-                Dpoints.Sort((x, y) => x.ParameterSet.Fitness.CompareTo(y.ParameterSet.Fitness));
+                Dpoints.Sort((x, y) => { int c = x.ParameterSet.Fitness.CompareTo(y.ParameterSet.Fitness); return c != 0 ? c : x.Index.CompareTo(y.Index); });
 
                 // Next reset the indexes, so that i = 0 represents the point With the smallest function value.
                 for (i = 0; i <= Dpoints.Count - 1; i++)
@@ -287,7 +287,7 @@ private void EvolveComplex(int alpha, ref List<PointFitness> Acomplex, double[]
                 {
 
                     // Rank the B points in terms of fitness (i.e., order of increasing function value.
-                    B.Sort((x, y) => x.ParameterSet.Fitness.CompareTo(y.ParameterSet.Fitness));
+                    B.Sort((x, y) => { int c = x.ParameterSet.Fitness.CompareTo(y.ParameterSet.Fitness); return c != 0 ? c : x.Index.CompareTo(y.Index); });
 
                     // Find centroid g excluding worst point and compute reflection of worst point about centroid r = 2g - u(worst)
                     Reflection(ref B, ref g, ref p);
@@ -350,7 +350,7 @@ private void EvolveComplex(int alpha, ref List<PointFitness> Acomplex, double[]
                     }
                 }
 
-                Acomplex.Sort((x, y) => x.ParameterSet.Fitness.CompareTo(y.ParameterSet.Fitness));
+                Acomplex.Sort((x, y) => { int c = x.ParameterSet.Fitness.CompareTo(y.ParameterSet.Fitness); return c != 0 ? c : x.Index.CompareTo(y.Index); });
             }
         }
 

diff --git a/Numerics/Sampling/MCMC/Base/MCMCSampler.cs b/Numerics/Sampling/MCMC/Base/MCMCSampler.cs
@@ -402,7 +402,7 @@ protected virtual ParameterSet[] InitializeChains()
                 return chainStates;
             }
 
-            var prng = new Random(PRNGSeed);
+            var prng = new MersenneTwister(PRNGSeed);
             var rnds = LatinHypercube.Random(InitialIterations, NumberOfParameters, prng.Next());
             var parameters = new double[NumberOfParameters];
             var tempPopulation = new List<ParameterSet>();       
@@ -634,14 +634,14 @@ public void Reset()
             _simulations = 0;
             MAPInitializationFailed = false;
             // Clear old memory and re-instantiate the result storage
-            _masterPRNG = new Random(PRNGSeed);
-            _chainPRNGs = new Random[NumberOfChains];
+            _masterPRNG = new MersenneTwister(PRNGSeed);
+            _chainPRNGs = new MersenneTwister[NumberOfChains];
             PopulationMatrix = new List<ParameterSet>();
             MarkovChains = new List<ParameterSet>[NumberOfChains];
             Output = new List<ParameterSet>[NumberOfChains];
             for (int i = 0; i < NumberOfChains; i++)
             {
-                _chainPRNGs[i] = new Random(_masterPRNG.Next());
+                _chainPRNGs[i] = new MersenneTwister(_masterPRNG.Next());
                 MarkovChains[i] = new List<ParameterSet>();
                 Output[i] = new List<ParameterSet>();
             }

diff --git a/Numerics/Utilities/ExtensionMethods.cs b/Numerics/Utilities/ExtensionMethods.cs
@@ -337,7 +337,7 @@ public static T[] RandomSubset<T>(this T[] array, int length, int seed = -1, boo
             if (length > array.Length)
                 throw new ArgumentException("The subset length must be less than or equal to the array length.");
 
-            Random random = seed > 0 ? new Random(seed) : new Random();
+            var random = seed > 0 ? new MersenneTwister(seed) : new MersenneTwister();
             var indexes = random.NextIntegers(0, array.Length, length, replace);
             var result = new T[indexes.Length];
             for (int i = 0; i < indexes.Length; i++)
@@ -484,7 +484,7 @@ public static void SetColumn<T>(this T[,] array, int index, T[] values)
             if (length > array.GetLength(0))
                 throw new ArgumentException("The subset length must be less than or equal to the array length.");
 
-            Random random = seed > 0 ? new Random(seed) : new Random();
+            var random = seed > 0 ? new MersenneTwister(seed) : new MersenneTwister();
             var indexes = random.NextIntegers(0, array.GetLength(0), length, replace);
             var result = new T[indexes.Length, array.GetLength(1)];
             for (int i = 0; i < indexes.Length; i++)
@@ -561,7 +561,7 @@ public static Vector RandomSubset(this Vector vector, int length, int seed = -1,
             if (length > vector.Length)
                 throw new ArgumentException("The subset length must be less than or equal to the vector length.");
 
-            Random random = seed > 0 ? new Random(seed) : new Random();
+            var random = seed > 0 ? new MersenneTwister(seed) : new MersenneTwister();
             var indexes = random.NextIntegers(0, vector.Length, length, replace);
             var result = new Vector(indexes.Length);
             for (int i = 0; i < indexes.Length; i++)
@@ -687,7 +687,7 @@ public static Matrix RandomSubset(this Matrix matrix, int length, int seed = -1,
             if (length > matrix.NumberOfRows)
                 throw new ArgumentException("The subset length must be less than or equal to the number of rows in the matrix.");
 
-            Random random = seed > 0 ? new Random(seed) : new Random();
+            var random = seed > 0 ? new MersenneTwister(seed) : new MersenneTwister();
             var indexes = random.NextIntegers(0, matrix.NumberOfRows, length, replace);
             var result = new Matrix(indexes.Length, matrix.NumberOfColumns);
             for (int i = 0; i < indexes.Length; i++)

diff --git a/Test_Numerics/AssemblyAttributes.cs b/Test_Numerics/AssemblyAttributes.cs
@@ -1,3 +1,3 @@
 using Microsoft.VisualStudio.TestTools.UnitTesting;
 
-[assembly: DoNotParallelize]
+[assembly: Parallelize(Scope = ExecutionScope.ClassLevel)]
diff --git a/Test_Numerics/Data/Statistics/Test_Probability.cs b/Test_Numerics/Data/Statistics/Test_Probability.cs
@@ -32,6 +32,7 @@
 using Numerics;
 using Numerics.Data.Statistics;
 using Numerics.Distributions;
+using Numerics.Sampling;
 using System;
 
 namespace Data.Statistics
@@ -107,7 +108,7 @@ public void Test_JointABCD_Independent_MVN()
               { r, r, 1, r },
               { r, r, r, 1 }};
 
-            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new Random(12345) };
+            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new MersenneTwister(12345) };
             var joint = Probability.JointProbabilityMVN(new double[] { A, B, C, D }, new int[] { 1, 1, 1, 1 }, mvn);
             Assert.AreEqual(0.021875, joint, 1E-6);
 
@@ -170,7 +171,7 @@ public void Test_JointABCD_PositivelyDependent_MVN()
               { r, r, 1, r },
               { r, r, r, 1 }};
 
-            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new Random(12345) };
+            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new MersenneTwister(12345) };
             var joint = Probability.JointProbabilityMVN(new double[] { A, B, C, D }, new int[] { 1, 1, 1, 1 }, mvn);
             Assert.AreEqual(0.25, joint, 1E-6);
 
@@ -233,7 +234,7 @@ public void Test_JointABCD_NegativelyDependent_MVN()
               { r, r, 1, r },
               { r, r, r, 1 }};
 
-            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new Random(12345) };
+            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new MersenneTwister(12345) };
             var joint = Probability.JointProbabilityMVN(new double[] { A, B, C, D }, new int[] { 1, 1, 1, 1 }, mvn);
             // True result comes from Monte Carlo simulation
             Assert.AreEqual(6.27974351606123E-14, joint, 1E-6);
@@ -306,7 +307,7 @@ public void Test_UnionABCD_Independent_MVN()
               { r, r, r, 1 }};
 
             // Get union from MVN
-            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new Random(12345) };
+            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new MersenneTwister(12345) };
             var union = Probability.UnionMVN(new double[] { A, B, C, D }, mvn);
             Assert.AreEqual(0.878125, union, 1E-6);
 
@@ -372,7 +373,7 @@ public void Test_UnionABCD_PositivelyDependent_MVN()
               { r, r, r, 1 }};
 
             // Get union from MVN
-            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new Random(12345) };
+            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new MersenneTwister(12345) };
             var union = Probability.UnionMVN(new double[] { A, B, C, D }, mvn);
             Assert.AreEqual(0.5, union, 1E-2);
         }
@@ -437,11 +438,11 @@ public void Test_UnionABCD_NegativelyDependent_MVN()
               { r, r, r, 1 }};
 
             // Get union from MVN
-            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new Random(12345) };
+            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new MersenneTwister(12345) };
             var union = Probability.UnionMVN(new double[] { A, B, C, D }, mvn);
 
             // True result comes from Monte Carlo simulation
-            Assert.AreEqual(0.985016583, union, 1E-5);
+            Assert.AreEqual(0.985016583, union, 1E-4);
         }
 
         /// <summary>
@@ -519,7 +520,7 @@ public void Test_ExclusiveABCD_Independent_MVN()
               { r, r, r, 1 }};
 
             // Get exclusive probabilities from MVN
-            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new Random(12345) };
+            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new MersenneTwister(12345) };
             var exclusive = Probability.ExclusiveMVN(new double[] { A, B, C, D }, mvn);
 
             // Test exclusive
@@ -612,7 +613,7 @@ public void Test_ExclusiveABCD_PositivelyDependent_MVN()
               { r, r, r, 1 }};
 
             // Get exclusive probabilities from MVN
-            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new Random(12345) };
+            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new MersenneTwister(12345) };
             var exclusive = Probability.ExclusiveMVN(new double[] { A, B, C, D }, mvn);
 
             // Both MVN and PCM have poor precision on this test
@@ -683,7 +684,7 @@ public void Test_ExclusiveABCD_NegativelyDependent_MVN()
               { r, r, r, 1 }};
 
             // Get exclusive probabilities from MVN
-            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new Random(12345) };
+            var mvn = new MultivariateNormal(new double[] { 0, 0, 0, 0 }, corr) { MVNUNI = new MersenneTwister(12345) };
             var exclusive = Probability.ExclusiveMVN(new double[] { A, B, C, D }, mvn);
 
             // Test exclusive

diff --git a/Test_Numerics/Data/Time Series/Test_TimeSeriesDownload.cs b/Test_Numerics/Data/Time Series/Test_TimeSeriesDownload.cs
@@ -50,6 +50,7 @@ namespace Data.TimeSeriesAnalysis
     /// </para>
     /// </remarks>
     [TestClass]
+    [DoNotParallelize]
     public class Test_TimeSeriesDownload
     {
         #region Station Lists

diff --git a/Test_Numerics/Distributions/Bivariate Copulas/Test_StudentTCopula.cs b/Test_Numerics/Distributions/Bivariate Copulas/Test_StudentTCopula.cs
@@ -33,6 +33,7 @@
 using Numerics.Data.Statistics;
 using Numerics.Distributions;
 using Numerics.Distributions.Copulas;
+using Numerics.Sampling;
 using Numerics.Sampling.MCMC;
 
 namespace Distributions.BivariateCopulas
@@ -238,7 +239,7 @@ public void Test_CDF_Monotonicity()
         public void Test_InverseCDF_RoundTrip()
         {
             var copula = new StudentTCopula(0.6, 5);
-            var rng = new Random(42);
+            var rng = new MersenneTwister(12345);
 
             for (int i = 0; i < 20; i++)
             {
@@ -261,7 +262,7 @@ public void Test_InverseCDF_Dependence()
             // With high positive ρ, InverseCDF(u, v) should produce correlated pairs
             var copula = new StudentTCopula(0.8, 5);
             int n = 1000;
-            var rng = new Random(42);
+            var rng = new MersenneTwister(12345);
             double sumProduct = 0;
             double sumU = 0, sumV = 0;
 

diff --git a/Test_Numerics/Distributions/Multivariate/Test_Multinomial.cs b/Test_Numerics/Distributions/Multivariate/Test_Multinomial.cs
@@ -31,6 +31,7 @@
 using System;
 using Microsoft.VisualStudio.TestTools.UnitTesting;
 using Numerics.Distributions;
+using Numerics.Sampling;
 
 namespace Distributions.Multivariate
 {
@@ -223,7 +224,7 @@ public void Test_Sampling()
         [TestMethod]
         public void Test_WeightedSample()
         {
-            var rng = new Random(42);
+            var rng = new MersenneTwister(12345);
             var weights = new double[] { 1.0, 3.0, 6.0 }; // 10%, 30%, 60%
             var counts = new int[3];
 
@@ -248,7 +249,7 @@ public void Test_WeightedSample()
         [TestMethod]
         public void Test_WeightedSample_EdgeCases()
         {
-            var rng = new Random(42);
+            var rng = new MersenneTwister(12345);
             // All weight on one category
             Assert.AreEqual(1, Multinomial.Sample(new double[] { 0.0, 1.0, 0.0 }, rng));
 

diff --git a/Test_Numerics/Mathematics/Optimization/Dynamic/BinaryHeapTesting.cs b/Test_Numerics/Mathematics/Optimization/Dynamic/BinaryHeapTesting.cs
@@ -33,6 +33,7 @@
 using System.Diagnostics;
 using Microsoft.VisualStudio.TestTools.UnitTesting;
 using Numerics.Mathematics.Optimization;
+using Numerics.Sampling;
 
 namespace Mathematics.Optimization
 {
@@ -72,7 +73,7 @@ public void HeapTest2()
             //Random Node weights
             float[] weights = new float[1000]; //= new float[] { .3f, .5f, 32f, 15f, 12f, .01f, -4f };
 
-            Random randy = new Random(42);
+            Random randy = new MersenneTwister(12345);
             for (int i = 0; i < weights.Length; i++)
             {
                 weights[i] = (float)randy.NextDouble();

diff --git a/Test_Numerics/Sampling/MCMC/Test_HMC.cs b/Test_Numerics/Sampling/MCMC/Test_HMC.cs
@@ -76,7 +76,7 @@ double logLH(double[] x)
 
             // Create and run MCMC sampler
             // NOTE: The HMC sampler is not adaptive, and the settings must be tuned by hand.
-            var sampler = new HMC(priors, logLH, stepSize: 2, steps: 10);
+            var sampler = new HMC(priors, logLH, stepSize: 2.5, steps: 10);
             sampler.Sample();
             var results = new MCMCResults(sampler);