USACE-RMC · HadenSmith · Jul 15, 2026 · Jul 10, 2026 · Jul 10, 2026 · Jul 11, 2026
diff --git a/.github/workflows/Snapshot.yml b/.github/workflows/Snapshot.yml
@@ -10,8 +10,8 @@ jobs:
     with:
       dotnet-version: '10.0.x'
       project-names: 'Numerics'
-      # Snapshot workflow appends .<run_number>-dev, so this tracks the next development patch after v2.1.2.
-      version: '2.1.3'
+      # Snapshot workflow appends .<run_number>-dev, so this tracks the next development patch after v2.1.3.
+      version: '2.1.4'
       # PR integration already runs the full multi-target test suite.
       run-tests: false
       nuget-source: 'https://www.hec.usace.army.mil/nexus/repository/consequences-nuget-public/'

diff --git a/.gitignore b/.gitignore
@@ -1,4 +1,4 @@
-## Ignore Visual Studio temporary files, build results, and    
+## Ignore Visual Studio temporary files, build results, and
 ## files generated by popular Visual Studio add-ons.
 ##
 ## Get latest from https://github.com/github/gitignore/blob/master/VisualStudio.gitignore
@@ -370,4 +370,11 @@ MigrationBackup/
 
 #Numerics Specific
 /TestResults
-/.claude/settings.local.json
+
+# Local developer guidance and settings. These files are for developer machines only.
+/CLAUDE.md
+/Claude.md
+/AGENTS.md
+/Agents.md
+/.claude/
+/.agents/
diff --git a/CITATION.cff b/CITATION.cff
@@ -2,8 +2,8 @@ cff-version: 1.2.0
 message: "If you use this software, please cite our article in the Journal of Open Source Software."
 type: software
 title: "Numerics: A .NET Library for Numerical Computing, Statistical Analysis, and Risk Assessment"
-version: "2.1.2"
-date-released: "2026-07-09"
+version: "2.1.3"
+date-released: "2026-07-15"
 license: 0BSD
 repository-code: "https://github.com/USACE-RMC/Numerics"
 url: "https://github.com/USACE-RMC/Numerics"

diff --git a/Directory.Build.props b/Directory.Build.props
@@ -0,0 +1,13 @@
+<Project>
+
+  <PropertyGroup>
+    <RepositoryRoot>$(MSBuildThisFileDirectory)</RepositoryRoot>
+    <EnforceXmlDocumentation Condition="'$(EnforceXmlDocumentation)' == ''">true</EnforceXmlDocumentation>
+  </PropertyGroup>
+
+  <PropertyGroup Condition="'$(EnforceXmlDocumentation)' == 'true'">
+    <GenerateDocumentationFile>true</GenerateDocumentationFile>
+    <WarningsAsErrors>$(WarningsAsErrors);CS1570;CS1571;CS1572;CS1573;CS1574;CS1584;CS1587;CS1589;CS1591</WarningsAsErrors>
+  </PropertyGroup>
+
+</Project>
diff --git a/Numerics/Data/Interpolation/Bilinear.cs b/Numerics/Data/Interpolation/Bilinear.cs
@@ -142,11 +142,11 @@ public double Interpolate(double x1, double x2)
             x1ub = X1Values[X1LI.Count - 1];
             if (X1Transform == Transform.Logarithmic)
             {
-                x1 = Math.Log10(x1);
-                x1i = Math.Log10(x1i);
-                x1ii = Math.Log10(x1ii);
-                x1lb = Math.Log10(x1lb);
-                x1ub = Math.Log10(x1ub);
+                x1 = Tools.Log10(x1);
+                x1i = Tools.Log10(x1i);
+                x1ii = Tools.Log10(x1ii);
+                x1lb = Tools.Log10(x1lb);
+                x1ub = Tools.Log10(x1ub);
             }
             else if (X1Transform == Transform.NormalZ)
             {
@@ -164,11 +164,11 @@ public double Interpolate(double x1, double x2)
             x2ub = X2Values[X2LI.Count - 1];
             if (X2Transform == Transform.Logarithmic)
             {
-                x2 = Math.Log10(x2);
-                x2j = Math.Log10(x2j);
-                x2jj = Math.Log10(x2jj);
-                x2lb = Math.Log10(x2lb);
-                x2ub = Math.Log10(x2ub);
+                x2 = Tools.Log10(x2);
+                x2j = Tools.Log10(x2j);
+                x2jj = Tools.Log10(x2jj);
+                x2lb = Tools.Log10(x2lb);
+                x2ub = Tools.Log10(x2ub);
             }
             else if (X2Transform == Transform.NormalZ)
             {
@@ -186,10 +186,10 @@ public double Interpolate(double x1, double x2)
             yijj = YValues[i, j + 1];
             if (YTransform == Transform.Logarithmic)
             {
-                yij = Math.Log10(yij);
-                yiij = Math.Log10(yiij);
-                yiijj = Math.Log10(yiijj);
-                yijj = Math.Log10(yijj);
+                yij = Tools.Log10(yij);
+                yiij = Tools.Log10(yiij);
+                yiijj = Tools.Log10(yiijj);
+                yijj = Tools.Log10(yijj);
             }
             else if (YTransform == Transform.NormalZ)
             {

diff --git a/Numerics/Data/Interpolation/Support/Search.cs b/Numerics/Data/Interpolation/Support/Search.cs
@@ -398,16 +398,34 @@ public static int Bisection(double x, IList<double> values, int start = 0, SortO
 
             // Perform bisection search
             int xlo = start, xhi = N, xm = 0;
-            while (xhi - xlo > 1)
+            if (order == SortOrder.Ascending)
             {
-                xm = xlo + (xhi - xlo >> 1); 
-                if (x >= values[xm] && order == SortOrder.Ascending)
+                while (xhi - xlo > 1)
                 {
-                    xlo = xm;
+                    xm = xlo + (xhi - xlo >> 1);
+                    if (x >= values[xm])
+                    {
+                        xlo = xm;
+                    }
+                    else
+                    {
+                        xhi = xm;
+                    }
                 }
-                else
+            }
+            else
+            {
+                while (xhi - xlo > 1)
                 {
-                    xhi = xm;
+                    xm = xlo + (xhi - xlo >> 1);
+                    if (x < values[xm])
+                    {
+                        xlo = xm;
+                    }
+                    else
+                    {
+                        xhi = xm;
+                    }
                 }
             }
             // Return XLO
@@ -480,16 +498,34 @@ public static int Bisection(double x, OrderedPairedData orderedPairedData, int s
 
             // Perform bisection search
             int xlo = start, xhi = N, xm = 0;
-            while (xhi - xlo > 1)
+            if (orderedPairedData.OrderX == SortOrder.Ascending)
             {
-                xm = xlo + (xhi - xlo >> 1);
-                if (x >= orderedPairedData[xm].X && orderedPairedData.OrderX == SortOrder.Ascending)
+                while (xhi - xlo > 1)
                 {
-                    xlo = xm;
+                    xm = xlo + (xhi - xlo >> 1);
+                    if (x >= orderedPairedData[xm].X)
+                    {
+                        xlo = xm;
+                    }
+                    else
+                    {
+                        xhi = xm;
+                    }
                 }
-                else
+            }
+            else
+            {
+                while (xhi - xlo > 1)
                 {
-                    xhi = xm;
+                    xm = xlo + (xhi - xlo >> 1);
+                    if (x < orderedPairedData[xm].X)
+                    {
+                        xlo = xm;
+                    }
+                    else
+                    {
+                        xhi = xm;
+                    }
                 }
             }
             // Return XLO
@@ -563,16 +599,34 @@ public static int Bisection(double x, IList<Ordinate> ordinates, int start = 0,
 
             // Perform bisection search
             int xlo = start, xhi = N, xm = 0;
-            while (xhi - xlo > 1)
+            if (order == SortOrder.Ascending)
             {
-                xm = xlo + (xhi - xlo >> 1);
-                if (x >= ordinates[xm].X && order == SortOrder.Ascending)
+                while (xhi - xlo > 1)
                 {
-                    xlo = xm;
+                    xm = xlo + (xhi - xlo >> 1);
+                    if (x >= ordinates[xm].X)
+                    {
+                        xlo = xm;
+                    }
+                    else
+                    {
+                        xhi = xm;
+                    }
                 }
-                else
+            }
+            else
+            {
+                while (xhi - xlo > 1)
                 {
-                    xhi = xm;
+                    xm = xlo + (xhi - xlo >> 1);
+                    if (x < ordinates[xm].X)
+                    {
+                        xlo = xm;
+                    }
+                    else
+                    {
+                        xhi = xm;
+                    }
                 }
             }
             // Return XLO
@@ -660,15 +714,15 @@ public static int Hunt(double x, IList<double> values, int start = 0, SortOrder
             }
             else
             {
-                if (x >= values[xlo] && order == SortOrder.Ascending)
+                if ((x >= values[xlo]) == (order == SortOrder.Ascending))
                 {
                     // Hunt up
                     for (;;)
                     {
                         // Not done hunting so double the increment
                         xhi = xlo + inc;
                         if (xhi >= N - 1) { xhi = N - 1; break; }
-                        else if (x < values[xhi] && order == SortOrder.Ascending) break;
+                        else if ((x < values[xhi]) == (order == SortOrder.Ascending)) break;
                         else
                         {
                             xlo = xhi;
@@ -684,7 +738,7 @@ public static int Hunt(double x, IList<double> values, int start = 0, SortOrder
                     {
                         xlo = xlo - inc;
                         if (xlo <= 0) { xlo = 0; break; }
-                        else if (x >= values[xlo] && order == SortOrder.Ascending) break;
+                        else if ((x >= values[xlo]) == (order == SortOrder.Ascending)) break;
                         else
                         {
                             xhi = xlo;
@@ -698,7 +752,7 @@ public static int Hunt(double x, IList<double> values, int start = 0, SortOrder
             while (xhi - xlo > 1)
             {
                 xm = xlo + (xhi - xlo >> 1);
-                if (x >= values[xm] && order == SortOrder.Ascending)
+                if ((x >= values[xm]) == (order == SortOrder.Ascending))
                 {
                     xlo = xm;
                 }
@@ -840,7 +894,7 @@ public static int Hunt(double xValue, OrderedPairedData orderedPairedData, int s
             while (XHI - XLO > 1)
             {
                 XM = XLO + (XHI - XLO >> 1);
-                if (X >= orderedPairedData[XM].X && ASCND == true)
+                if ((X >= orderedPairedData[XM].X) == ASCND)
                 {
                     XLO = XM;
                 }
@@ -988,7 +1042,7 @@ public static int Hunt(double xValue, IList<Ordinate> ordinateData, int start =
             while (XHI - XLO > 1)
             {
                 XM = XLO + (XHI - XLO >> 1);
-                if (X >= ordinateData[XM].X && ASCND == true)
+                if ((X >= ordinateData[XM].X) == ASCND)
                 {
                     XLO = XM;
                 }

diff --git a/Numerics/Data/Statistics/Histogram.cs b/Numerics/Data/Statistics/Histogram.cs
@@ -397,22 +397,32 @@ public void AddBin(Bin bin)
         public void AddData(double data)
         {
             SortBins();
-            int index = GetBinIndexOf(data);
             if (data <= LowerBound)
             {
-                _bins.First().LowerBound = data;
+                if (data < LowerBound)
+                {
+                    _bins.First().LowerBound = data;
+                    LowerBound = data;
+                    _areBinsSorted = false;
+                }
                 _bins.First().Frequency += 1;
-                _areBinsSorted = false;
             }
             else if (data >= UpperBound)
             {
-                _bins.Last().UpperBound = data;
+                if (data > UpperBound)
+                {
+                    _bins.Last().UpperBound = data;
+                    UpperBound = data;
+                }
                 _bins.Last().Frequency += 1;
-                _areBinsSorted = false;
             }
-            else if (index >= 0 && index < NumberOfBins)
+            else
             {
-                _bins[index].Frequency += 1;
+                int index = GetBinIndexOf(data);
+                if (index >= 0 && index < NumberOfBins)
+                {
+                    _bins[index].Frequency += 1;
+                }
             }
         }
 

diff --git a/Numerics/Data/Statistics/Probability.cs b/Numerics/Data/Statistics/Probability.cs
@@ -305,6 +305,7 @@ public static double JointProbabilityHPCM(IList<double> probabilities, int[] ind
 
             // Get z-values
             double zMin = -9, zMax = 9;
+            const double minimumCdf = 1E-300;
             var R = new double[n, n];
             Array.Copy(correlationMatrix, R, correlationMatrix.Length);
             int i, j, k, ir, ic;
@@ -325,7 +326,7 @@ public static double JointProbabilityHPCM(IList<double> probabilities, int[] ind
             z1 = R[0, 0];
             pdf = Normal.StandardPDF(z1);
             cdf = Normal.StandardCDF(z1);
-            //if (cdf < 1e-300) cdf = 1e-300;
+            if (cdf < minimumCdf) cdf = minimumCdf;
             A = pdf / cdf;
             B = A * (z1 + A);
             for (k = 1; k < n; k++)
@@ -352,7 +353,7 @@ public static double JointProbabilityHPCM(IList<double> probabilities, int[] ind
                 z1 = R[j, j - 1];
                 pdf = Normal.StandardPDF(z1);
                 cdf = Normal.StandardCDF(z1);
-                if (cdf < 1e-300) cdf = 1e-300;
+                if (cdf < minimumCdf) cdf = minimumCdf;
                 A = pdf / cdf;
                 B = A * (z1 + A);
                 for (k = j + 1; k < n; k++)

diff --git a/Numerics/Distributions/Bivariate Copulas/AMHCopula.cs b/Numerics/Distributions/Bivariate Copulas/AMHCopula.cs
@@ -82,6 +82,12 @@ public override double ThetaMaximum
         /// <inheritdoc/>
         public override ArgumentOutOfRangeException? ValidateParameter(double parameter, bool throwException)
         {
+            if (double.IsNaN(parameter) || double.IsInfinity(parameter))
+            {
+                var exception = new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter must be finite.");
+                if (throwException) throw exception;
+                return exception;
+            }
             if (parameter < ThetaMinimum)
             {
                 if (throwException) throw new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter θ (theta) must be greater than or equal to " + ThetaMinimum.ToString() + ".");
@@ -167,7 +173,7 @@ public override double[] InverseCDF(double u, double v)
         /// <inheritdoc/>
         public override BivariateCopula Clone()
         {
-            return new AMHCopula(Theta, MarginalDistributionX, MarginalDistributionY);
+            return new AMHCopula(Theta, CloneMarginal(MarginalDistributionX), CloneMarginal(MarginalDistributionY));
         }
 
         /// <summary>

diff --git a/Numerics/Distributions/Bivariate Copulas/Base/ArchimedeanCopula.cs b/Numerics/Distributions/Bivariate Copulas/Base/ArchimedeanCopula.cs
@@ -49,6 +49,12 @@ public override void SetCopulaParameters(double[] parameters)
         /// <inheritdoc/>
         public override ArgumentOutOfRangeException? ValidateParameter(double parameter, bool throwException)
         {
+            if (double.IsNaN(parameter) || double.IsInfinity(parameter))
+            {
+                var exception = new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter must be finite.");
+                if (throwException) throw exception;
+                return exception;
+            }
             if (parameter < ThetaMinimum)
             {
                 if (throwException) throw new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter θ (theta) must be greater than or equal to " + ThetaMinimum.ToString() + ".");
@@ -60,7 +66,11 @@ public override void SetCopulaParameters(double[] parameters)
                 if (throwException) throw new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter θ (theta) must be less than or equal to " + ThetaMaximum.ToString() + ".");
                 return new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter θ (theta) must be less than or equal to " + ThetaMaximum.ToString() + ".");
             }
-            return new ArgumentOutOfRangeException(nameof(Theta),"Parameter is valid");
+            // A valid parameter must return null: the Theta setter derives ParametersValid
+            // from 'ValidateParameter(value, false) is null', so returning a sentinel
+            // exception here left ParametersValid permanently false for every Archimedean
+            // family that did not override this method (Clayton, Gumbel, Joe).
+            return null;
         }
 
         /// <inheritdoc/>