diff --git a/.gitignore b/.gitignore
index f93f341..6c16898 100644
--- a/.gitignore
+++ b/.gitignore
@@ -32,3 +32,7 @@ _ReSharper*/
#Nuget packages folder
packages/
/.claude/settings.json
+
+# BenchmarkDotNet run artifacts (reports are copied into docs/benchmarks/ by hand)
+docs/benchmarks/artifacts/
+BenchmarkDotNet.Artifacts/
diff --git a/Numerics/NumericTest/FittingConditioningTests.cs b/Numerics/NumericTest/FittingConditioningTests.cs
new file mode 100644
index 0000000..0f82b53
--- /dev/null
+++ b/Numerics/NumericTest/FittingConditioningTests.cs
@@ -0,0 +1,195 @@
+using CSharpNumerics.Numerics.Objects;
+using CSharpNumerics.Statistics.Fitting;
+
+namespace NumericTest;
+
+///
+/// Accuracy tests for the least squares path on ill-conditioned designs.
+///
+///
+/// These exist to hold onto the reason the fitters were moved from the normal equations to QR.
+/// Forming XᵀX squares the condition number of X, and a Vandermonde design — any polynomial fit —
+/// is badly conditioned to begin with. The tolerances here are tight enough that a normal-equations
+/// solve cannot meet them, so if someone reverts the fitters to XᵀX these fail rather than
+/// silently losing digits.
+///
+[TestClass]
+public class FittingConditioningTests
+{
+ ///
+ /// Recovers the coefficients of a known polynomial sampled away from the origin, where the
+ /// Vandermonde matrix is badly conditioned.
+ ///
+ ///
+ /// Measured on this exact design: solving the normal equations recovers the coefficients to
+ /// a worst-case error of 5.6e-5, while QR on the design matrix reaches 7.2e-11 — close to six
+ /// decimal digits, which is what squaring a condition number costs. The 1e-9 tolerance sits
+ /// between the two, so this test fails if the fitters ever return to XᵀX.
+ ///
+ [TestMethod]
+ public void HighDegreePolynomial_RecoversKnownCoefficients()
+ {
+ // p(t) = 1 - 2t + 3t² - 4t³ + 5t⁴ - 6t⁵
+ var expected = new[] { 1.0, -2.0, 3.0, -4.0, 5.0, -6.0 };
+
+ var samples = 40;
+ var x = new double[samples];
+ var y = new double[samples];
+
+ for (var i = 0; i < samples; i++)
+ {
+ // Sampling over [1, 2] rather than [0, 1] already makes the columns strongly
+ // correlated; the normal equations square that.
+ var t = 1.0 + i / (double)(samples - 1);
+ x[i] = t;
+
+ var value = 0.0;
+ var power = 1.0;
+ for (var k = 0; k < expected.Length; k++)
+ {
+ value += expected[k] * power;
+ power *= t;
+ }
+
+ y[i] = value;
+ }
+
+ var result = LeastSquaresFitter.Fit(new VectorN(x), new VectorN(y), degree: 5);
+
+ for (var k = 0; k < expected.Length; k++)
+ {
+ Assert.AreEqual(expected[k], result.Coefficients[k], 1e-9,
+ $"Coefficient of t^{k}");
+ }
+ }
+
+ ///
+ /// A design matrix whose columns are nearly parallel. Squaring its condition number is
+ /// enough to lose the solution entirely; QR keeps it.
+ ///
+ [TestMethod]
+ public void NearlyCollinearColumns_StillRecoverTheSolution()
+ {
+ const double epsilon = 1e-7;
+ var samples = 30;
+
+ var design = new double[samples, 2];
+ var y = new double[samples];
+
+ for (var i = 0; i < samples; i++)
+ {
+ var t = i / (double)(samples - 1);
+
+ design[i, 0] = 1.0;
+ design[i, 1] = 1.0 + epsilon * t;
+
+ // Exact response for beta = (2, 3).
+ y[i] = 2.0 * design[i, 0] + 3.0 * design[i, 1];
+ }
+
+ var result = LeastSquaresFitter.Fit(design, new VectorN(y));
+
+ // Only the sum is well determined when the columns are this close, so that is what
+ // is asserted — but it must be right.
+ var fittedSum = result.Coefficients[0] + result.Coefficients[1];
+
+ Assert.AreEqual(5.0, fittedSum, 1e-6);
+ }
+
+ [TestMethod]
+ public void IllConditionedFit_ResidualsAreAtNoiseLevel()
+ {
+ var samples = 50;
+ var x = new double[samples];
+ var y = new double[samples];
+
+ for (var i = 0; i < samples; i++)
+ {
+ var t = 5.0 + 2.0 * i / (double)(samples - 1);
+ x[i] = t;
+ y[i] = 3.0 - 1.5 * t + 0.25 * t * t * t;
+ }
+
+ var result = LeastSquaresFitter.Fit(new VectorN(x), new VectorN(y), degree: 4);
+
+ // The model contains the truth, so the fit should be essentially exact.
+ for (var i = 0; i < samples; i++)
+ {
+ Assert.AreEqual(0.0, result.Residuals[i], 1e-8, $"Residual {i}");
+ }
+ }
+
+ ///
+ /// Standard errors now come from the triangular factor rather than from an explicitly
+ /// inverted Gram matrix. They must stay finite, positive and symmetric in the obvious way.
+ ///
+ [TestMethod]
+ public void StandardErrors_RemainWellDefinedOnAnIllConditionedDesign()
+ {
+ var samples = 40;
+ var x = new double[samples];
+ var y = new double[samples];
+ var random = new Random(99);
+
+ for (var i = 0; i < samples; i++)
+ {
+ var t = 8.0 + i / (double)(samples - 1);
+ x[i] = t;
+ y[i] = 1.0 + 0.5 * t - 0.1 * t * t + 0.01 * (random.NextDouble() - 0.5);
+ }
+
+ var result = LeastSquaresFitter.Fit(new VectorN(x), new VectorN(y), degree: 3);
+
+ for (var k = 0; k < result.StandardErrors.Length; k++)
+ {
+ Assert.IsFalse(double.IsNaN(result.StandardErrors[k]), $"SE {k} is NaN");
+ Assert.IsFalse(double.IsInfinity(result.StandardErrors[k]), $"SE {k} is infinite");
+ Assert.IsTrue(result.StandardErrors[k] >= 0.0, $"SE {k} is negative");
+ }
+ }
+
+ [TestMethod]
+ public void WeightedFit_WithWidelyDifferentWeights_RecoversTheSolution()
+ {
+ var samples = 25;
+ var design = new double[samples, 2];
+ var y = new double[samples];
+ var weights = new double[samples];
+
+ for (var i = 0; i < samples; i++)
+ {
+ var t = i / (double)(samples - 1);
+
+ design[i, 0] = 1.0;
+ design[i, 1] = t;
+ y[i] = 4.0 - 7.0 * t;
+
+ // Weights spanning eight orders of magnitude: squaring them in XᵀWX spans sixteen.
+ weights[i] = i % 2 == 0 ? 1e-4 : 1e4;
+ }
+
+ var result = WeightedLeastSquaresFitter.Fit(
+ design, new VectorN(y), new VectorN(weights));
+
+ Assert.AreEqual(4.0, result.Coefficients[0], 1e-8);
+ Assert.AreEqual(-7.0, result.Coefficients[1], 1e-8);
+ }
+
+ [TestMethod]
+ public void RankDeficientDesign_Throws()
+ {
+ // Two identical columns: no unique least squares solution exists.
+ var design = new double[5, 2];
+ var y = new double[5];
+
+ for (var i = 0; i < 5; i++)
+ {
+ design[i, 0] = i + 1;
+ design[i, 1] = i + 1;
+ y[i] = i;
+ }
+
+ Assert.ThrowsException(
+ () => LeastSquaresFitter.Fit(design, new VectorN(y)));
+ }
+}
diff --git a/Numerics/NumericTest/NaiveBayesTests.cs b/Numerics/NumericTest/NaiveBayesTests.cs
new file mode 100644
index 0000000..7425449
--- /dev/null
+++ b/Numerics/NumericTest/NaiveBayesTests.cs
@@ -0,0 +1,111 @@
+using CSharpNumerics.ML.Models.Classification;
+using CSharpNumerics.Numerics.Objects;
+
+namespace NumericTest;
+
+///
+/// Tests for the Gaussian naive Bayes classifier, which had no coverage before.
+///
+///
+/// Added alongside the fix to NumClasses, which threw
+/// NotImplementedException — the only such member left in the library.
+///
+[TestClass]
+public class NaiveBayesTests
+{
+ ///
+ /// Two well-separated Gaussian blobs: class 0 near the origin, class 1 far from it.
+ ///
+ private static (Matrix x, VectorN y) TwoBlobs()
+ {
+ var x = new double[,]
+ {
+ { 0.0, 0.0 }, { 0.2, -0.1 }, { -0.1, 0.2 }, { 0.1, 0.1 },
+ { 5.0, 5.0 }, { 5.2, 4.9 }, { 4.9, 5.2 }, { 5.1, 5.1 }
+ };
+
+ var y = new double[] { 0, 0, 0, 0, 1, 1, 1, 1 };
+
+ return (new Matrix(x), new VectorN(y));
+ }
+
+ [TestMethod]
+ public void NumClasses_IsSetByFit()
+ {
+ var (x, y) = TwoBlobs();
+
+ var model = new NaiveBayes();
+ model.Fit(x, y);
+
+ Assert.AreEqual(2, model.NumClasses);
+ }
+
+ [TestMethod]
+ public void NumClasses_CountsTheHighestLabelPlusOne()
+ {
+ // Labels 0 and 2 with nothing at 1: NumClasses sizes a label-indexed array,
+ // so it must be 3, matching the other classifiers.
+ var x = new Matrix(new double[,] { { 0.0 }, { 1.0 }, { 8.0 }, { 9.0 } });
+ var y = new VectorN(new double[] { 0, 0, 2, 2 });
+
+ var model = new NaiveBayes();
+ model.Fit(x, y);
+
+ Assert.AreEqual(3, model.NumClasses);
+ }
+
+ [TestMethod]
+ public void Predict_SeparatesTwoWellSeparatedBlobs()
+ {
+ var (x, y) = TwoBlobs();
+
+ var model = new NaiveBayes();
+ model.Fit(x, y);
+
+ var predictions = model.Predict(x);
+
+ for (var i = 0; i < y.Length; i++)
+ {
+ Assert.AreEqual(y[i], predictions[i], $"Sample {i}");
+ }
+ }
+
+ [TestMethod]
+ public void Predict_AssignsUnseenPointsToTheNearerClass()
+ {
+ var (x, y) = TwoBlobs();
+
+ var model = new NaiveBayes();
+ model.Fit(x, y);
+
+ var unseen = new Matrix(new double[,] { { 0.3, 0.3 }, { 4.7, 5.3 } });
+ var predictions = model.Predict(unseen);
+
+ Assert.AreEqual(0.0, predictions[0]);
+ Assert.AreEqual(1.0, predictions[1]);
+ }
+
+ [TestMethod]
+ public void PredictBeforeFit_Throws()
+ {
+ var model = new NaiveBayes();
+
+ Assert.ThrowsException(
+ () => model.Predict(new Matrix(new double[,] { { 1.0, 2.0 } })));
+ }
+
+ [TestMethod]
+ public void Clone_ReturnsAnUnfittedModel()
+ {
+ var (x, y) = TwoBlobs();
+
+ var model = new NaiveBayes();
+ model.Fit(x, y);
+
+ // Clone follows the convention of the other models: hyperparameters only, no fitted
+ // state, so grid search starts each fold from a clean estimator.
+ var clone = (NaiveBayes)model.Clone();
+
+ Assert.ThrowsException(() => clone.Predict(x));
+ }
+}
diff --git a/Numerics/NumericTest/NumericsTests.cs b/Numerics/NumericTest/NumericsTests.cs
index a0a3a8b..fa21752 100644
--- a/Numerics/NumericTest/NumericsTests.cs
+++ b/Numerics/NumericTest/NumericsTests.cs
@@ -219,7 +219,7 @@ public void FindRoots()
{
Func func = (double x) => Math.Pow(x,2) - 4;
var result = func.NewtonRaphson();
- Assert.IsTrue(Math.Abs(2) == 2);
+ Assert.AreEqual(2.0, result, 1e-6);
}
[TestMethod]
diff --git a/Numerics/NumericTest/PanelMethodTests.cs b/Numerics/NumericTest/PanelMethodTests.cs
new file mode 100644
index 0000000..cc2d5ae
--- /dev/null
+++ b/Numerics/NumericTest/PanelMethodTests.cs
@@ -0,0 +1,206 @@
+using CSharpNumerics.Physics.FluidDynamics.Aerodynamics;
+
+namespace NumericTest;
+
+///
+/// Characterisation tests for the 2-D panel method.
+///
+///
+/// Written to pin the solver's behaviour down before its hand-rolled Gaussian elimination is
+/// replaced with the shared LU decomposition.
+///
+/// These assert invariants that hold whatever the solver computes — symmetry, invariance of Cp
+/// under freestream scaling, net zero source strength, finiteness, sensitivity to incidence —
+/// rather than physical correctness. That is deliberate: the one test here that did check against
+/// a closed-form solution does not pass, and is left ignored with the question written down. So
+/// read this file as a change detector, not as validation that the physics is right.
+///
+[TestClass]
+public class PanelMethodTests
+{
+ ///
+ /// Builds a closed circular contour, traversed clockwise so the panel normals point
+ /// outward in the same sense the solver assumes for an airfoil.
+ ///
+ private static (double[] x, double[] y) Cylinder(int panels, double radius = 1.0)
+ {
+ var x = new double[panels + 1];
+ var y = new double[panels + 1];
+
+ for (var i = 0; i <= panels; i++)
+ {
+ var theta = -2.0 * Math.PI * i / panels;
+ x[i] = radius * Math.Cos(theta);
+ y[i] = radius * Math.Sin(theta);
+ }
+
+ return (x, y);
+ }
+
+ /// A symmetric diamond, the simplest closed body with sharp corners.
+ private static (double[] x, double[] y) Diamond()
+ {
+ var x = new double[] { 1.0, 0.5, 0.0, 0.5, 1.0 };
+ var y = new double[] { 0.0, -0.1, 0.0, 0.1, 0.0 };
+ return (x, y);
+ }
+
+ ///
+ /// Open question, not a passing assertion. Potential flow around a cylinder has
+ /// Cp = 1 − 4 sin²θ, ranging over [−3, 1]. This solver instead returns Cp ≈ 1 on every
+ /// panel (measured range [0.9999, 1.0000] at 200 panels, zero incidence), which means the
+ /// tangential velocity comes out ≈ 0 everywhere — not a potential flow around a body.
+ ///
+ ///
+ /// Two readings are possible and this test does not settle which: the solver may be wrong,
+ /// or a seamed cylinder may simply be outside its domain, since it imposes a Kutta condition
+ /// at a trailing edge that a cylinder does not have. Only the clockwise traversal was
+ /// measured. Left ignored rather than deleted so the question is not lost — investigating it
+ /// is separate work from the v4.3 migration.
+ ///
+ [TestMethod]
+ [Ignore("Unresolved: the solver returns Cp ~= 1 everywhere for a cylinder. See the summary above.")]
+ public void Cylinder_MatchesAnalyticPressureDistribution()
+ {
+ // Potential flow around a cylinder: Cp = 1 - 4 sin²θ, independent of radius.
+ var (x, y) = Cylinder(200);
+
+ var result = PanelMethod.Solve(x, y, alpha: 0.0);
+
+ for (var i = 0; i < result.Cp.Length; i++)
+ {
+ var theta = Math.Atan2(result.Ym[i], result.Xm[i]);
+ var expected = 1.0 - 4.0 * Math.Sin(theta) * Math.Sin(theta);
+
+ Assert.AreEqual(expected, result.Cp[i], 0.05,
+ $"Panel {i} at theta = {theta:F3}");
+ }
+ }
+
+ [TestMethod]
+ public void Cylinder_IsSymmetricAtZeroIncidence()
+ {
+ var (x, y) = Cylinder(120);
+
+ var result = PanelMethod.Solve(x, y, alpha: 0.0);
+
+ // A body symmetric about y = 0 at zero incidence must carry a symmetric
+ // pressure distribution.
+ for (var i = 0; i < result.Cp.Length; i++)
+ {
+ var mirror = NearestPanel(result, result.Xm[i], -result.Ym[i]);
+
+ Assert.AreEqual(result.Cp[i], result.Cp[mirror], 1e-6,
+ $"Panel {i} and its mirror {mirror}");
+ }
+ }
+
+ [TestMethod]
+ public void FreestreamScaling_LeavesPressureCoefficientUnchanged()
+ {
+ var (x, y) = Cylinder(80);
+
+ var unit = PanelMethod.Solve(x, y, alpha: 0.0, freestream: 1.0);
+ var scaled = PanelMethod.Solve(x, y, alpha: 0.0, freestream: 7.5);
+
+ // Cp is normalised by dynamic pressure, so it cannot depend on the freestream speed.
+ for (var i = 0; i < unit.Cp.Length; i++)
+ {
+ Assert.AreEqual(unit.Cp[i], scaled.Cp[i], 1e-9, $"Panel {i}");
+ }
+ }
+
+ [TestMethod]
+ public void SourceStrengths_SumToApproximatelyZeroForAClosedBody()
+ {
+ var (x, y) = Cylinder(160);
+
+ var result = PanelMethod.Solve(x, y, alpha: 0.0);
+
+ // A closed body in incompressible flow neither creates nor destroys mass.
+ var net = 0.0;
+ for (var i = 0; i < result.Sigma.Length; i++)
+ {
+ net += result.Sigma[i];
+ }
+
+ Assert.AreEqual(0.0, net, 1e-6);
+ }
+
+ [TestMethod]
+ public void Diamond_ProducesFiniteResultsOnASharpBody()
+ {
+ var (x, y) = Diamond();
+
+ var result = PanelMethod.Solve(x, y, alpha: 0.05);
+
+ Assert.AreEqual(4, result.Cp.Length);
+
+ foreach (var cp in result.Cp)
+ {
+ Assert.IsFalse(double.IsNaN(cp), "Cp must not be NaN");
+ Assert.IsFalse(double.IsInfinity(cp), "Cp must not be infinite");
+ }
+
+ foreach (var vt in result.Vt)
+ {
+ Assert.IsFalse(double.IsNaN(vt), "Vt must not be NaN");
+ }
+ }
+
+ [TestMethod]
+ public void IncidenceChangesThePressureDistribution()
+ {
+ var (x, y) = Cylinder(120);
+
+ var straight = PanelMethod.Solve(x, y, alpha: 0.0);
+ var inclined = PanelMethod.Solve(x, y, alpha: 0.2);
+
+ var changed = false;
+ for (var i = 0; i < straight.Cp.Length; i++)
+ {
+ if (Math.Abs(straight.Cp[i] - inclined.Cp[i]) > 1e-6)
+ {
+ changed = true;
+ break;
+ }
+ }
+
+ Assert.IsTrue(changed, "Angle of attack must affect the solution.");
+ }
+
+ [TestMethod]
+ public void TooFewPoints_Throws()
+ {
+ Assert.ThrowsException(
+ () => PanelMethod.Solve(new double[] { 0, 1, 0 }, new double[] { 0, 0, 0 }, 0.0));
+ }
+
+ [TestMethod]
+ public void NullCoordinates_Throws()
+ {
+ Assert.ThrowsException(
+ () => PanelMethod.Solve(null!, null!, 0.0));
+ }
+
+ private static int NearestPanel(PanelMethodResult result, double x, double y)
+ {
+ var best = 0;
+ var bestDistance = double.MaxValue;
+
+ for (var i = 0; i < result.Xm.Length; i++)
+ {
+ var dx = result.Xm[i] - x;
+ var dy = result.Ym[i] - y;
+ var distance = dx * dx + dy * dy;
+
+ if (distance < bestDistance)
+ {
+ bestDistance = distance;
+ best = i;
+ }
+ }
+
+ return best;
+ }
+}
diff --git a/Numerics/NumericTest/PcaTests.cs b/Numerics/NumericTest/PcaTests.cs
new file mode 100644
index 0000000..fcb6758
--- /dev/null
+++ b/Numerics/NumericTest/PcaTests.cs
@@ -0,0 +1,312 @@
+using CSharpNumerics.ML.DimensionalityReduction.Algorithms;
+using CSharpNumerics.Numerics.Objects;
+
+namespace NumericTest;
+
+///
+/// Tests for Principal Component Analysis.
+///
+///
+/// Written before the eigensolver behind PCA was swapped from power iteration with deflation to
+/// the shared EigenDecomposition. They assert mathematical properties of the
+/// decomposition — orthonormal components, variance ordering, alignment with known structure —
+/// so they hold for any correct eigensolver rather than pinning one implementation's output.
+///
+[TestClass]
+public class PcaTests
+{
+ ///
+ /// Points along the line y = 2x, so the leading principal direction is (1, 2) normalised.
+ ///
+ private static Matrix CollinearData()
+ {
+ var values = new double[7, 2];
+ for (var i = 0; i < 7; i++)
+ {
+ var t = i - 3.0;
+ values[i, 0] = t;
+ values[i, 1] = 2.0 * t;
+ }
+
+ return new Matrix(values);
+ }
+
+ /// Axis-aligned data whose variance is much larger along the first column.
+ private static Matrix AxisAlignedData()
+ {
+ var random = new Random(17);
+ var values = new double[60, 3];
+
+ for (var i = 0; i < 60; i++)
+ {
+ values[i, 0] = random.NextDouble() * 100.0;
+ values[i, 1] = random.NextDouble() * 10.0;
+ values[i, 2] = random.NextDouble();
+ }
+
+ return new Matrix(values);
+ }
+
+ [TestMethod]
+ public void Fit_CentersOnTheColumnMeans()
+ {
+ var data = new Matrix(new double[,] { { 1, 10 }, { 3, 20 }, { 5, 30 } });
+
+ var pca = new PCA { NComponents = 1, Seed = 1 };
+ pca.Fit(data);
+
+ Assert.AreEqual(3.0, pca.Mean[0], 1e-12);
+ Assert.AreEqual(20.0, pca.Mean[1], 1e-12);
+ }
+
+ [TestMethod]
+ public void FirstComponent_AlignsWithTheDirectionOfGreatestVariance()
+ {
+ var pca = new PCA { NComponents = 1, Seed = 1 };
+ pca.Fit(CollinearData());
+
+ // Expected direction (1, 2)/sqrt(5). Sign is arbitrary for an eigenvector, so compare
+ // the absolute value of the dot product with the unit expectation.
+ var expectedX = 1.0 / Math.Sqrt(5.0);
+ var expectedY = 2.0 / Math.Sqrt(5.0);
+
+ var alignment = Math.Abs(pca.Components[0, 0] * expectedX + pca.Components[0, 1] * expectedY);
+
+ Assert.AreEqual(1.0, alignment, 1e-6);
+ }
+
+ [TestMethod]
+ public void FirstComponent_PicksTheHighestVarianceAxis()
+ {
+ var pca = new PCA { NComponents = 1, Seed = 1 };
+ pca.Fit(AxisAlignedData());
+
+ // Variance is dominated by column 0, so the leading component must point essentially
+ // along it.
+ Assert.AreEqual(1.0, Math.Abs(pca.Components[0, 0]), 1e-3);
+ }
+
+ [TestMethod]
+ public void Components_AreUnitLength()
+ {
+ var pca = new PCA { NComponents = 3, Seed = 1 };
+ pca.Fit(AxisAlignedData());
+
+ for (var k = 0; k < pca.NComponents; k++)
+ {
+ var norm = 0.0;
+ for (var j = 0; j < 3; j++)
+ {
+ norm += pca.Components[k, j] * pca.Components[k, j];
+ }
+
+ Assert.AreEqual(1.0, Math.Sqrt(norm), 1e-8, $"Component {k}");
+ }
+ }
+
+ [TestMethod]
+ public void Components_AreMutuallyOrthogonal()
+ {
+ var pca = new PCA { NComponents = 3, Seed = 1 };
+ pca.Fit(AxisAlignedData());
+
+ for (var a = 0; a < pca.NComponents; a++)
+ {
+ for (var b = a + 1; b < pca.NComponents; b++)
+ {
+ var dot = 0.0;
+ for (var j = 0; j < 3; j++)
+ {
+ dot += pca.Components[a, j] * pca.Components[b, j];
+ }
+
+ Assert.AreEqual(0.0, dot, 1e-6, $"Components {a} and {b}");
+ }
+ }
+ }
+
+ [TestMethod]
+ public void ExplainedVariance_IsDescending()
+ {
+ var pca = new PCA { NComponents = 3, Seed = 1 };
+ pca.Fit(AxisAlignedData());
+
+ for (var k = 1; k < pca.NComponents; k++)
+ {
+ Assert.IsTrue(
+ pca.ExplainedVariance[k] <= pca.ExplainedVariance[k - 1] + 1e-9,
+ $"Component {k} explains more variance than {k - 1}: " +
+ $"{pca.ExplainedVariance[k]} > {pca.ExplainedVariance[k - 1]}");
+ }
+ }
+
+ [TestMethod]
+ public void ExplainedVarianceRatio_SumsToOneWhenAllComponentsAreKept()
+ {
+ var pca = new PCA { NComponents = 3, Seed = 1 };
+ pca.Fit(AxisAlignedData());
+
+ var total = 0.0;
+ for (var k = 0; k < pca.NComponents; k++)
+ {
+ total += pca.ExplainedVarianceRatio[k];
+ }
+
+ Assert.AreEqual(1.0, total, 1e-6);
+ }
+
+ [TestMethod]
+ public void RankDeficientData_PutsAllVarianceInTheFirstComponent()
+ {
+ var pca = new PCA { NComponents = 2, Seed = 1 };
+ pca.Fit(CollinearData());
+
+ // The data lies on a line, so the second component explains nothing.
+ Assert.AreEqual(1.0, pca.ExplainedVarianceRatio[0], 1e-6);
+ Assert.AreEqual(0.0, pca.ExplainedVarianceRatio[1], 1e-6);
+ }
+
+ [TestMethod]
+ public void Transform_ProjectsOntoTheRequestedNumberOfComponents()
+ {
+ var data = AxisAlignedData();
+
+ var pca = new PCA { NComponents = 2, Seed = 1 };
+ var reduced = pca.FitTransform(data);
+
+ Assert.AreEqual(data.rowLength, reduced.rowLength);
+ Assert.AreEqual(2, reduced.columnLength);
+ }
+
+ [TestMethod]
+ public void Transform_ProducesUncorrelatedZeroMeanScores()
+ {
+ var data = AxisAlignedData();
+
+ var pca = new PCA { NComponents = 3, Seed = 1 };
+ var scores = pca.FitTransform(data);
+
+ // Projections onto principal directions are centred and uncorrelated by construction.
+ for (var k = 0; k < 3; k++)
+ {
+ var mean = 0.0;
+ for (var i = 0; i < scores.rowLength; i++)
+ {
+ mean += scores.values[i, k];
+ }
+
+ Assert.AreEqual(0.0, mean / scores.rowLength, 1e-8, $"Score column {k} mean");
+ }
+
+ for (var a = 0; a < 3; a++)
+ {
+ for (var b = a + 1; b < 3; b++)
+ {
+ var covariance = 0.0;
+ for (var i = 0; i < scores.rowLength; i++)
+ {
+ covariance += scores.values[i, a] * scores.values[i, b];
+ }
+
+ Assert.AreEqual(0.0, covariance / scores.rowLength, 1e-6,
+ $"Score columns {a} and {b}");
+ }
+ }
+ }
+
+ [TestMethod]
+ public void Transform_ReconstructsFullRankDataWhenAllComponentsAreKept()
+ {
+ var data = AxisAlignedData();
+
+ var pca = new PCA { NComponents = 3, Seed = 1 };
+ var scores = pca.FitTransform(data);
+
+ // X ≈ mean + scores · components
+ for (var i = 0; i < data.rowLength; i++)
+ {
+ for (var j = 0; j < 3; j++)
+ {
+ var reconstructed = pca.Mean[j];
+ for (var k = 0; k < 3; k++)
+ {
+ reconstructed += scores.values[i, k] * pca.Components[k, j];
+ }
+
+ Assert.AreEqual(data.values[i, j], reconstructed, 1e-8, $"Element ({i}, {j})");
+ }
+ }
+ }
+
+ [TestMethod]
+ public void WideData_UsesTheDualPathAndStillProducesUnitComponents()
+ {
+ // Fewer samples than features drives the Gram-matrix branch.
+ var random = new Random(5);
+ var values = new double[4, 12];
+ for (var i = 0; i < 4; i++)
+ {
+ for (var j = 0; j < 12; j++)
+ {
+ values[i, j] = random.NextDouble();
+ }
+ }
+
+ var pca = new PCA { NComponents = 2, Seed = 1 };
+ pca.Fit(new Matrix(values));
+
+ Assert.AreEqual(2, pca.NComponents);
+
+ for (var k = 0; k < 2; k++)
+ {
+ var norm = 0.0;
+ for (var j = 0; j < 12; j++)
+ {
+ norm += pca.Components[k, j] * pca.Components[k, j];
+ }
+
+ Assert.AreEqual(1.0, Math.Sqrt(norm), 1e-8, $"Component {k}");
+ }
+ }
+
+ [TestMethod]
+ public void NComponents_IsClampedToTheAvailableRank()
+ {
+ var data = new Matrix(new double[,] { { 1, 2 }, { 3, 5 } });
+
+ var pca = new PCA { NComponents = 10, Seed = 1 };
+ pca.Fit(data);
+
+ Assert.AreEqual(2, pca.NComponents);
+ }
+
+ [TestMethod]
+ public void TransformBeforeFit_Throws()
+ {
+ var pca = new PCA { NComponents = 1 };
+
+ Assert.ThrowsException(
+ () => pca.Transform(AxisAlignedData()));
+ }
+
+ [TestMethod]
+ public void ZeroComponents_Throws()
+ {
+ var pca = new PCA { NComponents = 0 };
+
+ Assert.ThrowsException(() => pca.Fit(AxisAlignedData()));
+ }
+
+ [TestMethod]
+ public void Clone_CarriesTheFittedState()
+ {
+ var pca = new PCA { NComponents = 2, Seed = 1 };
+ pca.Fit(AxisAlignedData());
+
+ var clone = (PCA)pca.Clone();
+
+ Assert.AreEqual(pca.NComponents, clone.NComponents);
+ Assert.AreEqual(pca.Components[0, 0], clone.Components[0, 0], 1e-15);
+ Assert.AreEqual(pca.Mean[0], clone.Mean[0], 1e-15);
+ }
+}
diff --git a/Numerics/NumericTest/RootFindingTests.cs b/Numerics/NumericTest/RootFindingTests.cs
new file mode 100644
index 0000000..0afc697
--- /dev/null
+++ b/Numerics/NumericTest/RootFindingTests.cs
@@ -0,0 +1,346 @@
+using CSharpNumerics.Numerics;
+using CSharpNumerics.Numerics.RootFinding;
+
+namespace NumericTest;
+
+[TestClass]
+public class RootFindingTests
+{
+ // f(x) = x² − 4, root at x = 2.
+ private static readonly Func Quadratic = x => x * x - 4.0;
+
+ // f(x) = cos(x) − x, root at the Dottie number.
+ private static readonly Func CosineMinusX = x => Math.Cos(x) - x;
+ private const double DottieNumber = 0.7390851332151607;
+
+ // f(x) = x³ − 2x − 5, the classic test case from Numerical Recipes.
+ private static readonly Func Cubic = x => x * x * x - 2.0 * x - 5.0;
+ private const double CubicRoot = 2.0945514815423265;
+
+ #region Bisection
+
+ [TestMethod]
+ public void Bisection_Quadratic_FindsRoot()
+ {
+ var result = RootFinder.Bisection(Quadratic, 0.0, 5.0);
+
+ Assert.IsTrue(result.Converged, result.ToString());
+ Assert.AreEqual(2.0, result.Value, 1e-9);
+ }
+
+ [TestMethod]
+ public void Bisection_ReversedBracket_IsAccepted()
+ {
+ var result = RootFinder.Bisection(Quadratic, 5.0, 0.0);
+
+ Assert.IsTrue(result.Converged);
+ Assert.AreEqual(2.0, result.Value, 1e-9);
+ }
+
+ [TestMethod]
+ public void Bisection_RootExactlyOnEndpoint_ReturnsImmediately()
+ {
+ var result = RootFinder.Bisection(Quadratic, 2.0, 5.0);
+
+ Assert.IsTrue(result.Converged);
+ Assert.AreEqual(2.0, result.Value, 0.0);
+ Assert.AreEqual(0, result.Iterations);
+ }
+
+ [TestMethod]
+ public void Bisection_NoSignChange_Throws()
+ {
+ Assert.ThrowsException(
+ () => RootFinder.Bisection(Quadratic, 3.0, 5.0));
+ }
+
+ [TestMethod]
+ public void Bisection_NonSmoothFunction_StillConverges()
+ {
+ // |x| − 1 has a corner at the origin but a clean sign change at x = 1.
+ Func kinked = x => Math.Abs(x) - 1.0;
+
+ var result = RootFinder.Bisection(kinked, 0.5, 3.0);
+
+ Assert.IsTrue(result.Converged);
+ Assert.AreEqual(1.0, result.Value, 1e-9);
+ }
+
+ #endregion
+
+ #region Secant
+
+ [TestMethod]
+ public void Secant_Quadratic_FindsRoot()
+ {
+ var result = RootFinder.Secant(Quadratic, 0.0, 5.0);
+
+ Assert.IsTrue(result.Converged, result.ToString());
+ Assert.AreEqual(2.0, result.Value, 1e-9);
+ }
+
+ [TestMethod]
+ public void Secant_EqualStartingPoints_Throws()
+ {
+ Assert.ThrowsException(
+ () => RootFinder.Secant(Quadratic, 1.0, 1.0));
+ }
+
+ [TestMethod]
+ public void Secant_HorizontalSecant_ReportsFailureInsteadOfDividingByZero()
+ {
+ // A constant non-zero function gives f(x₁) − f(x₀) = 0: no step is defined.
+ Func constant = _ => 3.0;
+
+ var result = RootFinder.Secant(constant, 0.0, 1.0);
+
+ Assert.IsFalse(result.Converged);
+ Assert.IsFalse(double.IsNaN(result.Value));
+ }
+
+ #endregion
+
+ #region Brent
+
+ [TestMethod]
+ public void Brent_Quadratic_FindsRoot()
+ {
+ var result = RootFinder.Brent(Quadratic, 0.0, 5.0);
+
+ Assert.IsTrue(result.Converged, result.ToString());
+ Assert.AreEqual(2.0, result.Value, 1e-12);
+ }
+
+ [TestMethod]
+ public void Brent_Transcendental_MatchesKnownRoot()
+ {
+ var result = RootFinder.Brent(CosineMinusX, 0.0, 1.0);
+
+ Assert.IsTrue(result.Converged, result.ToString());
+ Assert.AreEqual(DottieNumber, result.Value, 1e-12);
+ }
+
+ [TestMethod]
+ public void Brent_Cubic_MatchesKnownRoot()
+ {
+ var result = RootFinder.Brent(Cubic, 2.0, 3.0);
+
+ Assert.IsTrue(result.Converged, result.ToString());
+ Assert.AreEqual(CubicRoot, result.Value, 1e-12);
+ }
+
+ [TestMethod]
+ public void Brent_TripleRoot_ConvergesDespiteVanishingDerivative()
+ {
+ // (x − 1)³ is flat at its root, which defeats plain interpolation.
+ Func tripleRoot = x => Math.Pow(x - 1.0, 3);
+
+ var result = RootFinder.Brent(tripleRoot, -2.0, 4.0);
+
+ Assert.IsTrue(result.Converged, result.ToString());
+ Assert.AreEqual(1.0, result.Value, 1e-4);
+ }
+
+ [TestMethod]
+ public void Brent_CubeRoot_ConvergesWhereNewtonDiverges()
+ {
+ // ∛x has an infinite derivative at the root — Newton's method never converges on it.
+ Func cubeRoot = Math.Cbrt;
+
+ var brent = RootFinder.Brent(cubeRoot, -1.0, 2.0);
+ var newton = RootFinder.Newton(cubeRoot, 1.0);
+
+ Assert.IsTrue(brent.Converged, brent.ToString());
+ Assert.AreEqual(0.0, brent.Value, 1e-9);
+ Assert.IsFalse(newton.Converged, "Newton is expected to fail on the cube root.");
+ }
+
+ [TestMethod]
+ public void Brent_NoSignChange_Throws()
+ {
+ Assert.ThrowsException(
+ () => RootFinder.Brent(Quadratic, 3.0, 5.0));
+ }
+
+ [TestMethod]
+ public void Brent_ConvergesInFarFewerIterationsThanBisection()
+ {
+ var brent = RootFinder.Brent(CosineMinusX, 0.0, 1.0);
+ var bisection = RootFinder.Bisection(CosineMinusX, 0.0, 1.0);
+
+ Assert.IsTrue(brent.Converged);
+ Assert.IsTrue(bisection.Converged);
+ Assert.IsTrue(
+ brent.Iterations < bisection.Iterations / 2,
+ $"Brent took {brent.Iterations} iterations, bisection took {bisection.Iterations}.");
+ }
+
+ #endregion
+
+ #region Newton
+
+ [TestMethod]
+ public void Newton_Quadratic_FindsRoot()
+ {
+ var result = RootFinder.Newton(Quadratic, 1.0);
+
+ Assert.IsTrue(result.Converged, result.ToString());
+ Assert.AreEqual(2.0, result.Value, 1e-6);
+ }
+
+ [TestMethod]
+ public void Newton_AnalyticDerivative_IsMoreAccurateThanFiniteDifference()
+ {
+ var analytic = RootFinder.Newton(Quadratic, x => 2.0 * x, 1.0);
+ var finiteDifference = RootFinder.Newton(Quadratic, 1.0);
+
+ Assert.IsTrue(analytic.Converged, analytic.ToString());
+ Assert.IsTrue(finiteDifference.Converged, finiteDifference.ToString());
+ Assert.AreEqual(2.0, analytic.Value, 1e-12);
+ Assert.IsTrue(
+ analytic.Residual <= finiteDifference.Residual,
+ $"Analytic residual {analytic.Residual:G6} should not exceed " +
+ $"finite-difference residual {finiteDifference.Residual:G6}.");
+ }
+
+ [TestMethod]
+ public void Newton_StopsEarlyInsteadOfRunningTheFullBudget()
+ {
+ var result = RootFinder.Newton(Quadratic, x => 2.0 * x, 1.0, maxIterations: 100);
+
+ Assert.IsTrue(result.Converged);
+ Assert.IsTrue(result.Iterations < 20, $"Took {result.Iterations} iterations.");
+ }
+
+ [TestMethod]
+ public void Newton_VanishingDerivative_ReportsFailureInsteadOfNaN()
+ {
+ // x³ − 3x + 3 has f'(1) = 0 exactly, so no Newton step exists from x₀ = 1.
+ Func function = x => x * x * x - 3.0 * x + 3.0;
+ Func derivative = x => 3.0 * x * x - 3.0;
+
+ var result = RootFinder.Newton(function, derivative, 1.0);
+
+ Assert.IsFalse(result.Converged);
+ Assert.IsFalse(double.IsNaN(result.Value));
+ Assert.IsFalse(double.IsInfinity(result.Value));
+ }
+
+ [TestMethod]
+ public void Newton_DivergentStart_ReportsFailure()
+ {
+ // arctan diverges under Newton's method for any |x₀| greater than ≈1.3917.
+ var result = RootFinder.Newton(Math.Atan, x => 1.0 / (1.0 + x * x), 5.0);
+
+ Assert.IsFalse(result.Converged);
+ }
+
+ [TestMethod]
+ public void Newton_NoRealRoot_ReportsFailure()
+ {
+ Func noRoot = x => x * x + 1.0;
+
+ var result = RootFinder.Newton(noRoot, x => 2.0 * x, 3.0);
+
+ Assert.IsFalse(result.Converged);
+ }
+
+ [TestMethod]
+ public void Newton_NullDerivative_Throws()
+ {
+ Assert.ThrowsException(
+ () => RootFinder.Newton(Quadratic, null!, 1.0));
+ }
+
+ #endregion
+
+ #region Argument validation
+
+ [TestMethod]
+ public void NonPositiveTolerance_Throws()
+ {
+ Assert.ThrowsException(
+ () => RootFinder.Brent(Quadratic, 0.0, 5.0, tolerance: 0.0));
+ }
+
+ [TestMethod]
+ public void ZeroIterationBudget_Throws()
+ {
+ Assert.ThrowsException(
+ () => RootFinder.Brent(Quadratic, 0.0, 5.0, maxIterations: 0));
+ }
+
+ [TestMethod]
+ public void DegenerateBracket_Throws()
+ {
+ Assert.ThrowsException(
+ () => RootFinder.Brent(Quadratic, 2.5, 2.5));
+ }
+
+ [TestMethod]
+ public void NullFunction_Throws()
+ {
+ Assert.ThrowsException(
+ () => RootFinder.Brent(null!, 0.0, 5.0));
+ }
+
+ #endregion
+
+ #region RootResult
+
+ [TestMethod]
+ public void EnsureConverged_OnSuccess_ReturnsValue()
+ {
+ var result = RootFinder.Brent(Quadratic, 0.0, 5.0);
+
+ Assert.AreEqual(2.0, result.EnsureConverged(), 1e-12);
+ }
+
+ [TestMethod]
+ public void EnsureConverged_OnFailure_Throws()
+ {
+ var result = RootFinder.Newton(Math.Atan, x => 1.0 / (1.0 + x * x), 5.0);
+
+ Assert.ThrowsException(() => result.EnsureConverged());
+ }
+
+ [TestMethod]
+ public void Result_ReportsResidualAtReturnedPoint()
+ {
+ var result = RootFinder.Brent(Cubic, 2.0, 3.0);
+
+ Assert.AreEqual(Math.Abs(Cubic(result.Value)), result.Residual, 1e-15);
+ }
+
+ #endregion
+
+ #region FindRoot facade and NewtonRaphson compatibility
+
+ [TestMethod]
+ public void FindRoot_DelegatesToBrent()
+ {
+ var result = CosineMinusX.FindRoot(0.0, 1.0);
+
+ Assert.IsTrue(result.Converged);
+ Assert.AreEqual(DottieNumber, result.Value, 1e-12);
+ Assert.AreEqual("Brent", result.Method);
+ }
+
+ [TestMethod]
+ public void NewtonRaphson_KeepsWorkingThroughTheNewImplementation()
+ {
+ var root = Quadratic.NewtonRaphson();
+
+ Assert.AreEqual(2.0, root, 1e-6);
+ }
+
+ [TestMethod]
+ public void NewtonRaphson_RespectsInitialGuess()
+ {
+ var negativeRoot = Quadratic.NewtonRaphson(-1.0);
+
+ Assert.AreEqual(-2.0, negativeRoot, 1e-6);
+ }
+
+ #endregion
+}
diff --git a/Numerics/Numerics.Benchmarks/DenseLinearAlgebraBenchmarks.cs b/Numerics/Numerics.Benchmarks/DenseLinearAlgebraBenchmarks.cs
new file mode 100644
index 0000000..35ff2ba
--- /dev/null
+++ b/Numerics/Numerics.Benchmarks/DenseLinearAlgebraBenchmarks.cs
@@ -0,0 +1,118 @@
+using BenchmarkDotNet.Attributes;
+using CSharpNumerics.Numerics.LinearAlgebra;
+using CSharpNumerics.Numerics.LinearAlgebra.Decompositions;
+using CSharpNumerics.Numerics.Objects;
+
+namespace Numerics.Benchmarks;
+
+///
+/// Baseline for the dense linear algebra kernels: the triple-loop multiply and the
+/// decompositions built in v4.1.
+///
+///
+/// These are the hot paths that v4.3 routes more callers through and that v4.4 plans to
+/// vectorise, so this is the measurement both are judged against.
+///
+[MemoryDiagnoser]
+public class DenseLinearAlgebraBenchmarks
+{
+ [Params(64, 256)]
+ public int N;
+
+ private Matrix general = default!;
+ private Matrix other = default!;
+ private Matrix symmetricPositiveDefinite = default!;
+ private VectorN rightHandSide = default!;
+ private LuDecomposition luFactorization = default!;
+
+ [GlobalSetup]
+ public void Setup()
+ {
+ // Fixed seed: the same matrices every run, so results are comparable across commits.
+ var random = new Random(20260930);
+
+ general = RandomMatrix(random, N, N);
+ other = RandomMatrix(random, N, N);
+ symmetricPositiveDefinite = MakeSymmetricPositiveDefinite(random, N);
+
+ var values = new double[N];
+ for (var i = 0; i < N; i++)
+ {
+ values[i] = random.NextDouble();
+ }
+
+ rightHandSide = new VectorN(values);
+ luFactorization = general.Lu();
+ }
+
+ [Benchmark(Description = "Matrix × Matrix")]
+ public Matrix MatrixMultiply() => general * other;
+
+ [Benchmark(Description = "Matrix × Vector")]
+ public VectorN MatrixVectorMultiply() => general * rightHandSide;
+
+ [Benchmark(Description = "LU factorize")]
+ public LuDecomposition LuFactorize() => general.Lu();
+
+ [Benchmark(Description = "QR factorize")]
+ public QrDecomposition QrFactorize() => general.Qr();
+
+ [Benchmark(Description = "Cholesky factorize")]
+ public CholeskyDecomposition CholeskyFactorize() => symmetricPositiveDefinite.Cholesky();
+
+ /// Factor and solve — what a caller pays who does not keep the factorization.
+ [Benchmark(Description = "LU factorize + solve")]
+ public VectorN LuFactorizeAndSolve() => general.Lu().Solve(rightHandSide);
+
+ ///
+ /// Solve only, against a factorization computed once. The gap between this and
+ /// is what the v4.3 migration is meant to recover
+ /// wherever several right-hand sides share one matrix.
+ ///
+ [Benchmark(Description = "LU solve (reusing factorization)")]
+ public VectorN LuSolveReusingFactorization() => luFactorization.Solve(rightHandSide);
+
+ [Benchmark(Description = "Matrix.Inverse")]
+ public Matrix Inverse() => general.Inverse();
+
+ private static Matrix RandomMatrix(Random random, int rows, int columns)
+ {
+ var values = new double[rows, columns];
+ for (var i = 0; i < rows; i++)
+ {
+ for (var j = 0; j < columns; j++)
+ {
+ values[i, j] = random.NextDouble() * 2.0 - 1.0;
+ }
+ }
+
+ return new Matrix(values);
+ }
+
+ ///
+ /// Builds BᵀB + nI, which is symmetric and diagonally dominant, so Cholesky is defined.
+ ///
+ private static Matrix MakeSymmetricPositiveDefinite(Random random, int n)
+ {
+ var b = RandomMatrix(random, n, n);
+ var values = new double[n, n];
+
+ for (var i = 0; i < n; i++)
+ {
+ for (var j = 0; j < n; j++)
+ {
+ var sum = 0.0;
+ for (var k = 0; k < n; k++)
+ {
+ sum += b.values[k, i] * b.values[k, j];
+ }
+
+ values[i, j] = sum;
+ }
+
+ values[i, i] += n;
+ }
+
+ return new Matrix(values);
+ }
+}
diff --git a/Numerics/Numerics.Benchmarks/MachineLearningBenchmarks.cs b/Numerics/Numerics.Benchmarks/MachineLearningBenchmarks.cs
new file mode 100644
index 0000000..9304328
--- /dev/null
+++ b/Numerics/Numerics.Benchmarks/MachineLearningBenchmarks.cs
@@ -0,0 +1,72 @@
+using BenchmarkDotNet.Attributes;
+using CSharpNumerics.ML.Models.Regression;
+using CSharpNumerics.Numerics.Objects;
+
+namespace Numerics.Benchmarks;
+
+///
+/// Baseline for the neural network training loop — the reference point for the
+/// allocation-free training work planned for v4.4.
+///
+///
+/// One epoch over a fixed synthetic regression set. The allocation column matters more
+/// here than the timing: every forward and backward pass currently allocates new
+/// and instances.
+///
+[MemoryDiagnoser]
+public class MachineLearningBenchmarks
+{
+ [Params(256, 2048)]
+ public int Samples;
+
+ private const int Features = 16;
+
+ private Matrix inputs = default!;
+ private VectorN targets = default!;
+
+ [GlobalSetup]
+ public void Setup()
+ {
+ var random = new Random(20260930);
+
+ var x = new double[Samples, Features];
+ var y = new double[Samples];
+
+ for (var i = 0; i < Samples; i++)
+ {
+ var sum = 0.0;
+ for (var j = 0; j < Features; j++)
+ {
+ var value = random.NextDouble() * 2.0 - 1.0;
+ x[i, j] = value;
+ sum += value * (j + 1);
+ }
+
+ y[i] = sum / Features;
+ }
+
+ inputs = new Matrix(x);
+ targets = new VectorN(y);
+ }
+
+ ///
+ /// A single training epoch of a 16 → 32 → 16 → 1 network. Early stopping is held off
+ /// so the measurement is one full pass rather than a variable number of them.
+ ///
+ [Benchmark(Description = "MLP regressor, one training epoch")]
+ public MLPRegressor TrainOneEpoch()
+ {
+ var model = new MLPRegressor
+ {
+ HiddenLayers = new[] { 32, 16 },
+ Epochs = 1,
+ BatchSize = 32,
+ LearningRate = 0.01,
+ Patience = int.MaxValue
+ };
+
+ model.Fit(inputs, targets);
+
+ return model;
+ }
+}
diff --git a/Numerics/Numerics.Benchmarks/Numerics.Benchmarks.csproj b/Numerics/Numerics.Benchmarks/Numerics.Benchmarks.csproj
new file mode 100644
index 0000000..b7cd8ff
--- /dev/null
+++ b/Numerics/Numerics.Benchmarks/Numerics.Benchmarks.csproj
@@ -0,0 +1,27 @@
+
+
+
+ Exe
+ net10.0
+ latest
+ enable
+ enable
+
+
+ false
+
+
+ Release
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Numerics/Numerics.Benchmarks/Program.cs b/Numerics/Numerics.Benchmarks/Program.cs
new file mode 100644
index 0000000..d577af4
--- /dev/null
+++ b/Numerics/Numerics.Benchmarks/Program.cs
@@ -0,0 +1,12 @@
+using System.Reflection;
+using BenchmarkDotNet.Running;
+
+// Entry point for the benchmark suite.
+//
+// dotnet run -c Release --project Numerics/Numerics.Benchmarks (menu)
+// dotnet run -c Release --project Numerics/Numerics.Benchmarks -- --filter *Dense*
+// dotnet run -c Release --project Numerics/Numerics.Benchmarks -- --list flat
+//
+// Results are written to BenchmarkDotNet.Artifacts/; copy the markdown report into
+// docs/benchmarks/ when recording a baseline.
+BenchmarkSwitcher.FromAssembly(Assembly.GetExecutingAssembly()).Run(args);
diff --git a/Numerics/Numerics.Benchmarks/README.md b/Numerics/Numerics.Benchmarks/README.md
new file mode 100644
index 0000000..8e57504
--- /dev/null
+++ b/Numerics/Numerics.Benchmarks/README.md
@@ -0,0 +1,65 @@
+# Numerics.Benchmarks
+
+Measurement project for CSharpNumerics, built on [BenchmarkDotNet](https://benchmarkdotnet.org/).
+
+It exists so that performance claims can be checked rather than asserted. The guiding
+principle in the roadmap is *measure before optimising*: the SIMD and allocation-free work
+planned for v4.4 is judged against the baseline recorded here, and the v4.3 migration of
+call sites onto the shared decompositions is expected to leave these numbers no worse.
+
+## Why it is a separate project
+
+The core library has no external dependencies, and it stays that way. BenchmarkDotNet lives
+here and nowhere else. The project is `IsPackable=false`, so it never reaches the NuGet
+package, and CI runs tests against `NumericTest` directly, so benchmarks never run there —
+they are far too slow for per-commit execution.
+
+## Running
+
+Release configuration is required; BenchmarkDotNet refuses to measure a Debug build.
+
+```powershell
+# Pick from a menu
+dotnet run -c Release --project Numerics/Numerics.Benchmarks
+
+# Everything (the full baseline — takes a while)
+dotnet run -c Release --project Numerics/Numerics.Benchmarks -- --filter *
+
+# One group
+dotnet run -c Release --project Numerics/Numerics.Benchmarks -- --filter *Dense*
+dotnet run -c Release --project Numerics/Numerics.Benchmarks -- --filter *Sparse*
+dotnet run -c Release --project Numerics/Numerics.Benchmarks -- --filter *MachineLearning*
+
+# List without running
+dotnet run -c Release --project Numerics/Numerics.Benchmarks -- --list flat
+
+# Smoke test that every benchmark executes (one cold iteration each, timings meaningless)
+dotnet run -c Release --project Numerics/Numerics.Benchmarks -- --job Dry --filter *
+```
+
+Reports land in `docs/benchmarks/artifacts/` (git-ignored). To record a baseline, copy the
+generated markdown report into `docs/benchmarks/` under a name that says what it measured and
+at which commit.
+
+## What is measured
+
+| Group | Benchmarks | Why |
+|-------|-----------|-----|
+| `DenseLinearAlgebraBenchmarks` | Matrix×Matrix, Matrix×Vector, LU / QR / Cholesky factorization, LU factor-and-solve, LU solve reusing a factorization, `Matrix.Inverse` | The triple-loop multiply is the hot path in neural network training, FEM and every decomposition. The two LU solve variants show what reusing a factorization is worth. |
+| `SparseLinearAlgebraBenchmarks` | SpMV, PCG solve | The 5-point Laplacian on an interior grid — the shape a 2-D FEM or finite difference problem assembles. SpMV is the inner loop of every iterative solver. |
+| `MachineLearningBenchmarks` | One MLP training epoch | Reference point for the allocation-free training loops planned for v4.4. |
+
+Every benchmark carries `[MemoryDiagnoser]`. For this library the **Allocated** column is often
+more informative than the timing: the current kernels allocate a fresh `Matrix` or `VectorN`
+for every intermediate result, and that is precisely what later work aims to remove.
+
+## Writing a benchmark
+
+- Seed any randomness with a fixed constant in `[GlobalSetup]` so runs are comparable across
+ commits.
+- Return the result from the benchmark method. A method that returns `void` and computes into
+ a local can be optimised away entirely.
+- Keep setup out of the measured method unless the setup *is* what you are measuring — the
+ two LU benchmarks are a deliberate pair that separates factoring from solving.
+- Build the problem correctly. A conjugate gradient benchmark on a matrix that is accidentally
+ non-symmetric measures a solver that has no reason to converge, not the solver's real cost.
diff --git a/Numerics/Numerics.Benchmarks/SparseLinearAlgebraBenchmarks.cs b/Numerics/Numerics.Benchmarks/SparseLinearAlgebraBenchmarks.cs
new file mode 100644
index 0000000..bae1f62
--- /dev/null
+++ b/Numerics/Numerics.Benchmarks/SparseLinearAlgebraBenchmarks.cs
@@ -0,0 +1,99 @@
+using BenchmarkDotNet.Attributes;
+using CSharpNumerics.Numerics.Objects;
+
+namespace Numerics.Benchmarks;
+
+///
+/// Baseline for the sparse path: SpMV and the preconditioned conjugate gradient solver
+/// that Assembler2D runs on.
+///
+///
+/// The system is the 5-point Laplacian on a × mesh —
+/// the same shape a 2-D finite element or finite difference problem produces, so the
+/// numbers transfer to real FEM workloads.
+///
+[MemoryDiagnoser]
+public class SparseLinearAlgebraBenchmarks
+{
+ [Params(32, 128)]
+ public int Grid;
+
+ private SparseMatrix laplacian = default!;
+ private VectorN rightHandSide = default!;
+
+ /// Degrees of freedom in the assembled system — the interior nodes.
+ public int Unknowns => (Grid - 2) * (Grid - 2);
+
+ [GlobalSetup]
+ public void Setup()
+ {
+ laplacian = BuildLaplacian(Grid);
+
+ var values = new double[Unknowns];
+ for (var i = 0; i < values.Length; i++)
+ {
+ values[i] = 1.0;
+ }
+
+ rightHandSide = new VectorN(values);
+ }
+
+ [Benchmark(Description = "SpMV (sparse matrix × vector)")]
+ public VectorN SparseMatrixVectorMultiply() => laplacian.Multiply(rightHandSide);
+
+ ///
+ /// Conjugate gradient with a Jacobi preconditioner. Iteration count grows with the
+ /// grid, so this scales worse than SpMV alone — that growth is what stronger
+ /// preconditioners would attack.
+ ///
+ [Benchmark(Description = "PCG solve")]
+ public VectorN SolvePcg() => laplacian.SolvePCG(rightHandSide, tolerance: 1e-8);
+
+ ///
+ /// Assembles the 5-point Laplacian over the interior nodes only, with the Dirichlet
+ /// boundary folded into the right-hand side rather than kept as rows.
+ ///
+ ///
+ /// Eliminating the boundary rather than pinning it with identity rows is what keeps the
+ /// matrix symmetric: an identity row zeroes a row but leaves the matching column entry
+ /// in its interior neighbours, and conjugate gradient has no convergence guarantee on a
+ /// non-symmetric system.
+ ///
+ private static SparseMatrix BuildLaplacian(int grid)
+ {
+ var interior = grid - 2;
+ var triplets = new List<(int row, int col, double val)>();
+ int Index(int x, int y) => y * interior + x;
+
+ for (var y = 0; y < interior; y++)
+ {
+ for (var x = 0; x < interior; x++)
+ {
+ var row = Index(x, y);
+ triplets.Add((row, row, 4.0));
+
+ if (x > 0)
+ {
+ triplets.Add((row, Index(x - 1, y), -1.0));
+ }
+
+ if (x < interior - 1)
+ {
+ triplets.Add((row, Index(x + 1, y), -1.0));
+ }
+
+ if (y > 0)
+ {
+ triplets.Add((row, Index(x, y - 1), -1.0));
+ }
+
+ if (y < interior - 1)
+ {
+ triplets.Add((row, Index(x, y + 1), -1.0));
+ }
+ }
+ }
+
+ return SparseMatrix.FromTriplets(interior * interior, interior * interior, triplets);
+ }
+}
diff --git a/Numerics/Numerics.sln b/Numerics/Numerics.sln
index 8722a23..c2bd1fb 100644
--- a/Numerics/Numerics.sln
+++ b/Numerics/Numerics.sln
@@ -9,6 +9,8 @@ Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "NumericTest", "NumericTest\
EndProject
Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Numerics", "Numerics", "{35409048-A951-1400-ED33-5FB42D2B3E23}"
EndProject
+Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Numerics.Benchmarks", "Numerics.Benchmarks\Numerics.Benchmarks.csproj", "{2CD34815-7653-4FEC-9069-210E4DC8723B}"
+EndProject
Global
GlobalSection(SolutionConfigurationPlatforms) = preSolution
Debug|Any CPU = Debug|Any CPU
@@ -43,6 +45,18 @@ Global
{26D48AFB-281B-452F-B8C8-5BCCD105925C}.Release|x64.Build.0 = Release|Any CPU
{26D48AFB-281B-452F-B8C8-5BCCD105925C}.Release|x86.ActiveCfg = Release|Any CPU
{26D48AFB-281B-452F-B8C8-5BCCD105925C}.Release|x86.Build.0 = Release|Any CPU
+ {2CD34815-7653-4FEC-9069-210E4DC8723B}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+ {2CD34815-7653-4FEC-9069-210E4DC8723B}.Debug|Any CPU.Build.0 = Debug|Any CPU
+ {2CD34815-7653-4FEC-9069-210E4DC8723B}.Debug|x64.ActiveCfg = Debug|Any CPU
+ {2CD34815-7653-4FEC-9069-210E4DC8723B}.Debug|x64.Build.0 = Debug|Any CPU
+ {2CD34815-7653-4FEC-9069-210E4DC8723B}.Debug|x86.ActiveCfg = Debug|Any CPU
+ {2CD34815-7653-4FEC-9069-210E4DC8723B}.Debug|x86.Build.0 = Debug|Any CPU
+ {2CD34815-7653-4FEC-9069-210E4DC8723B}.Release|Any CPU.ActiveCfg = Release|Any CPU
+ {2CD34815-7653-4FEC-9069-210E4DC8723B}.Release|Any CPU.Build.0 = Release|Any CPU
+ {2CD34815-7653-4FEC-9069-210E4DC8723B}.Release|x64.ActiveCfg = Release|Any CPU
+ {2CD34815-7653-4FEC-9069-210E4DC8723B}.Release|x64.Build.0 = Release|Any CPU
+ {2CD34815-7653-4FEC-9069-210E4DC8723B}.Release|x86.ActiveCfg = Release|Any CPU
+ {2CD34815-7653-4FEC-9069-210E4DC8723B}.Release|x86.Build.0 = Release|Any CPU
EndGlobalSection
GlobalSection(SolutionProperties) = preSolution
HideSolutionNode = FALSE
diff --git a/Numerics/Numerics/CSharpNumerics.csproj b/Numerics/Numerics/CSharpNumerics.csproj
index fc03a3c..1edde34 100644
--- a/Numerics/Numerics/CSharpNumerics.csproj
+++ b/Numerics/Numerics/CSharpNumerics.csproj
@@ -1,4 +1,4 @@
-
+
net10.0;net8.0;netstandard2.1
latest
@@ -20,10 +20,6 @@
-
-
-
-
diff --git a/Numerics/Numerics/ML/DimensionalityReduction/Algorithms/PCA.cs b/Numerics/Numerics/ML/DimensionalityReduction/Algorithms/PCA.cs
index 6b54332..aa76fd3 100644
--- a/Numerics/Numerics/ML/DimensionalityReduction/Algorithms/PCA.cs
+++ b/Numerics/Numerics/ML/DimensionalityReduction/Algorithms/PCA.cs
@@ -1,5 +1,6 @@
-using CSharpNumerics.ML.DimensionalityReduction.Interfaces;
+using CSharpNumerics.ML.DimensionalityReduction.Interfaces;
using CSharpNumerics.ML.Models.Interfaces;
+using CSharpNumerics.Numerics.LinearAlgebra.Decompositions;
using CSharpNumerics.Numerics.Objects;
using System;
using System.Collections.Generic;
@@ -8,15 +9,30 @@ namespace CSharpNumerics.ML.DimensionalityReduction.Algorithms;
///
/// Principal Component Analysis (PCA).
-/// Projects data onto the top eigenvectors of the covariance matrix.
-/// Uses the power iteration method with deflation for eigendecomposition.
+/// Projects data onto the top eigenvectors of the covariance matrix,
+/// found with .
///
public class PCA : IDimensionalityReducer, IHasHyperparameters
{
// ── Hyperparameters ──────────────────────────────────────────
public int NComponents { get; set; } = 2;
+
+ ///
+ /// No longer affects the result. Retained so existing calls and grid-search parameter
+ /// dictionaries keep working.
+ ///
+ ///
+ /// PCA used to find eigenpairs by power iteration with deflation, which needed an iteration
+ /// budget, a convergence tolerance and a random starting vector. It now uses
+ /// , which is direct and deterministic, so none of the three
+ /// has anything left to control.
+ ///
public int MaxIterations { get; set; } = 1000;
+
+ ///
public double Tolerance { get; set; } = 1e-8;
+
+ ///
public int? Seed { get; set; }
// ── Results (available after Fit) ────────────────────────────
@@ -67,58 +83,49 @@ public void Fit(Matrix X)
Components = new double[effectiveComponents, d];
ExplainedVariance = new double[effectiveComponents];
- var rng = Seed.HasValue ? new Random(Seed.Value) : new Random();
-
if (n >= d)
{
- // Standard PCA: covariance matrix (d × d)
- double[,] cov = ComputeCovariance(centered, n, d);
+ // Standard PCA: eigendecomposition of the covariance matrix (d × d).
+ var eigen = new EigenDecomposition(new Matrix(ComputeCovariance(centered, n, d)));
+ var eigenvalues = eigen.RealEigenvalues;
+ var eigenvectors = eigen.EigenVectors;
+ // Eigenvalues come back ascending, so the leading components are at the end.
for (int k = 0; k < effectiveComponents; k++)
{
- var (eigenvalue, eigenvector) = PowerIteration(cov, d, rng);
- ExplainedVariance[k] = eigenvalue;
+ int index = d - 1 - k;
+ ExplainedVariance[k] = eigenvalues[index];
for (int j = 0; j < d; j++)
- Components[k, j] = eigenvector[j];
-
- // Deflate
- for (int i = 0; i < d; i++)
- for (int j = 0; j < d; j++)
- cov[i, j] -= eigenvalue * eigenvector[i] * eigenvector[j];
+ Components[k, j] = eigenvectors.values[j, index];
}
}
else
{
- // Dual PCA: Gram matrix (n × n) — used when n << d to avoid OOM
- double[,] gram = ComputeGramMatrix(centered, n, d);
+ // Dual PCA: the Gram matrix (n × n) when n < d, which shares the covariance
+ // matrix's non-zero eigenvalues — both are already scaled by 1/(n-1).
+ var eigen = new EigenDecomposition(new Matrix(ComputeGramMatrix(centered, n, d)));
+ var eigenvalues = eigen.RealEigenvalues;
+ var eigenvectors = eigen.EigenVectors;
for (int k = 0; k < effectiveComponents; k++)
{
- var (eigenvalue, eigenvector) = PowerIteration(gram, n, rng);
+ int index = n - 1 - k;
- // Convert eigenvalue: gram eigenvalue = (n-1) * covariance eigenvalue
- double covEigenvalue = eigenvalue;
-
- // Recover d-dimensional eigenvector: v = X^T * u, then normalize
+ // Recover the d-dimensional direction: v = Xᵀu, then normalise.
var component = new double[d];
for (int j = 0; j < d; j++)
{
double sum = 0;
for (int i = 0; i < n; i++)
- sum += centered[i, j] * eigenvector[i];
+ sum += centered[i, j] * eigenvectors.values[i, index];
component[j] = sum;
}
NormalizeArray(component);
- ExplainedVariance[k] = covEigenvalue;
+ ExplainedVariance[k] = eigenvalues[index];
for (int j = 0; j < d; j++)
Components[k, j] = component[j];
-
- // Deflate gram matrix
- for (int i = 0; i < n; i++)
- for (int j = 0; j < n; j++)
- gram[i, j] -= eigenvalue * eigenvector[i] * eigenvector[j];
}
}
@@ -264,51 +271,6 @@ private static double ComputeTotalVariance(double[,] centered, int n, int d)
return total;
}
- private (double eigenvalue, double[] eigenvector) PowerIteration(
- double[,] matrix, int d, Random rng)
- {
- // Random initial vector
- var v = new double[d];
- for (int i = 0; i < d; i++)
- v[i] = rng.NextDouble() - 0.5;
- Normalize(v);
-
- double eigenvalue = 0;
-
- for (int iter = 0; iter < MaxIterations; iter++)
- {
- // w = A * v
- var w = new double[d];
- for (int i = 0; i < d; i++)
- {
- double sum = 0;
- for (int j = 0; j < d; j++)
- sum += matrix[i, j] * v[j];
- w[i] = sum;
- }
-
- // Rayleigh quotient: eigenvalue = v^T * w
- double newEigenvalue = 0;
- for (int i = 0; i < d; i++)
- newEigenvalue += v[i] * w[i];
-
- // Normalize
- Normalize(w);
-
- if (Math.Abs(newEigenvalue - eigenvalue) < Tolerance)
- {
- eigenvalue = newEigenvalue;
- v = w;
- break;
- }
-
- eigenvalue = newEigenvalue;
- v = w;
- }
-
- return (eigenvalue, v);
- }
-
private static void Normalize(double[] v)
{
double norm = 0;
diff --git a/Numerics/Numerics/ML/Models/Classification/NaiveBayes.cs b/Numerics/Numerics/ML/Models/Classification/NaiveBayes.cs
index 005950c..0954c09 100644
--- a/Numerics/Numerics/ML/Models/Classification/NaiveBayes.cs
+++ b/Numerics/Numerics/ML/Models/Classification/NaiveBayes.cs
@@ -15,7 +15,12 @@ public class NaiveBayes : IClassificationModel
private int _numFeatures;
private bool _fitted;
- public int NumClasses => throw new NotImplementedException();
+ ///
+ /// Number of classes, taken as the highest label seen during plus one,
+ /// matching the other implementations so callers can
+ /// use it to size label-indexed arrays.
+ ///
+ public int NumClasses { get; private set; }
public void Fit(Matrix X, VectorN y)
{
@@ -25,6 +30,7 @@ public void Fit(Matrix X, VectorN y)
_priors = new();
var classes = y.Values.Distinct().Select(v => (int)v).ToArray();
+ NumClasses = (int)y.Values.Max() + 1;
foreach (var cls in classes)
{
diff --git a/Numerics/Numerics/ML/ReinforcementLearning/Algorithms/Tabular/SARSA.cs b/Numerics/Numerics/ML/ReinforcementLearning/Algorithms/Tabular/SARSA.cs
index dcba4e3..5b99081 100644
--- a/Numerics/Numerics/ML/ReinforcementLearning/Algorithms/Tabular/SARSA.cs
+++ b/Numerics/Numerics/ML/ReinforcementLearning/Algorithms/Tabular/SARSA.cs
@@ -1,4 +1,4 @@
-using CSharpNumerics.ML.ReinforcementLearning.Core;
+using CSharpNumerics.ML.ReinforcementLearning.Core;
using CSharpNumerics.ML.ReinforcementLearning.Interfaces;
using CSharpNumerics.Numerics.Objects;
using System;
@@ -18,7 +18,6 @@ public class SARSA : TabularAgent
public override string Name => "SARSA";
private Transition _pending;
- private int _pendingNextAction;
private bool _hasPending;
public SARSA(int numStates, int numActions, Func stateMapper)
diff --git a/Numerics/Numerics/Numerics/DifferentialEquationExtensions.cs b/Numerics/Numerics/Numerics/DifferentialEquationExtensions.cs
index a07ef38..657429a 100644
--- a/Numerics/Numerics/Numerics/DifferentialEquationExtensions.cs
+++ b/Numerics/Numerics/Numerics/DifferentialEquationExtensions.cs
@@ -581,70 +581,20 @@ public static Vector GaussElimination(this Matrix matrix, Vector vector)
///
/// Solves a linear system using Gaussian elimination with partial pivoting.
///
- /// Coefficient matrix (will be modified in-place).
- /// Right-hand side vector (will be modified in-place).
+ /// Coefficient matrix A.
+ /// Right-hand side vector b.
/// Solution vector as a list.
+ ///
+ /// Delegates to , which performs the same elimination with the
+ /// same partial pivoting. Two behaviours changed with that move, both for the better:
+ /// the arguments are no longer modified in place, and a solution component that is exactly
+ /// zero now comes back as zero. The previous back substitution returned the pivot
+ /// A[i, i] in that case, so any system with a zero in its solution was answered
+ /// wrongly without complaint.
+ ///
public static List GaussElimination(this Matrix matrix, List vector)
{
- var n = vector.Count;
-
- for (var p = 0; p < n; p++)
- {
- var max = p;
- for (var i = p + 1; i < n; i++)
- {
- if (Math.Abs(matrix.values[i, p]) > Math.Abs(matrix.values[max, p]))
- {
- max = i;
- }
- }
-
- for (var c = 0; c < matrix.columnLength; c++)
- {
- var temp = matrix.values[p, c];
-
- matrix.values[p, c] = matrix.values[max, c];
-
- matrix.values[max, c] = temp;
- }
-
- var t = vector[p];
- vector[p] = vector[max];
- vector[max] = t;
-
- for (var i = p + 1; i < n; i++)
- {
- var alpha = matrix.values[i, p] / matrix.values[p, p];
- vector[i] -= alpha * vector[p];
- for (var j = p; j < n; j++)
- {
- matrix.values[i, j] -= alpha * matrix.values[p, j];
- }
- }
- }
-
- var x = new double[n];
- for (var i = n - 1; i >= 0; i--)
- {
- var sum = 0.0;
- for (var j = i + 1; j < n; j++)
- {
- sum += matrix.values[i, j] * x[j];
- }
- if (matrix.values[i, i] != 0)
- {
- if ((vector[i] - sum) != 0)
- {
- x[i] = (vector[i] - sum) / matrix.values[i, i];
- }
-
- if ((vector[i] - sum) == 0)
- {
- x[i] = matrix.values[i, i];
- }
- }
- }
- return x.ToList();
+ return new LuDecomposition(matrix).Solve(vector);
}
///
diff --git a/Numerics/Numerics/Numerics/FiniteElement/Assembler1D.cs b/Numerics/Numerics/Numerics/FiniteElement/Assembler1D.cs
index adfdf89..c2beca7 100644
--- a/Numerics/Numerics/Numerics/FiniteElement/Assembler1D.cs
+++ b/Numerics/Numerics/Numerics/FiniteElement/Assembler1D.cs
@@ -1,8 +1,9 @@
-namespace CSharpNumerics.Numerics.FiniteElement;
+namespace CSharpNumerics.Numerics.FiniteElement;
using System;
using System.Collections.Generic;
using CSharpNumerics.Numerics.FiniteElement.Interfaces;
+using CSharpNumerics.Numerics.LinearAlgebra.Decompositions;
using CSharpNumerics.Numerics.Objects;
///
@@ -76,7 +77,7 @@ public void ApplyNodalLoad(int nodeIndex, int localDof, double value)
///
/// Solves the global system Ku = F after applying Dirichlet boundary conditions
- /// by row/column elimination and Gaussian elimination with partial pivoting.
+ /// by row/column elimination, then solving densely via LU with partial pivoting.
///
/// Dictionary mapping global DOF index → prescribed value.
/// Solution vector of all DOFs.
@@ -106,45 +107,7 @@ public VectorN Solve(Dictionary fixedDofs)
b[dof] = val;
}
- // Gaussian elimination with partial pivoting
- for (int col = 0; col < n; col++)
- {
- int maxRow = col;
- double maxVal = Math.Abs(A[col, col]);
- for (int row = col + 1; row < n; row++)
- {
- double v = Math.Abs(A[row, col]);
- if (v > maxVal) { maxVal = v; maxRow = row; }
- }
-
- if (maxRow != col)
- {
- for (int j = col; j < n; j++)
- (A[col, j], A[maxRow, j]) = (A[maxRow, j], A[col, j]);
- (b[col], b[maxRow]) = (b[maxRow], b[col]);
- }
-
- double pivot = A[col, col];
- for (int row = col + 1; row < n; row++)
- {
- double factor = A[row, col] / pivot;
- for (int j = col; j < n; j++)
- A[row, j] -= factor * A[col, j];
- b[row] -= factor * b[col];
- }
- }
-
- // Back substitution
- var x = new double[n];
- for (int i = n - 1; i >= 0; i--)
- {
- double sum = b[i];
- for (int j = i + 1; j < n; j++)
- sum -= A[i, j] * x[j];
- x[i] = sum / A[i, i];
- }
-
- return new VectorN(x);
+ return new LuDecomposition(new Matrix(A)).Solve(new VectorN(b));
}
private int[] GetElementDofs(int nodeA, int nodeB)
diff --git a/Numerics/Numerics/Numerics/Interpolation/CubicSplineInterpolation.cs b/Numerics/Numerics/Numerics/Interpolation/CubicSplineInterpolation.cs
index b22bf8f..0ae748e 100644
--- a/Numerics/Numerics/Numerics/Interpolation/CubicSplineInterpolation.cs
+++ b/Numerics/Numerics/Numerics/Interpolation/CubicSplineInterpolation.cs
@@ -1,5 +1,7 @@
-using System;
+using System;
using System.Linq;
+using CSharpNumerics.Numerics.LinearAlgebra.Decompositions;
+using CSharpNumerics.Numerics.Objects;
namespace CSharpNumerics.Numerics.Interpolation;
@@ -244,38 +246,13 @@ private void ComputeNotAKnot(double[] h, double[] M)
A[m, m - 2] = h[m - 1]; A[m, m - 1] = -(h[m - 2] + h[m - 1]); A[m, m] = h[m - 2];
b[m] = 0;
- // Gaussian elimination with partial pivoting
- for (int col = 0; col < n; col++)
- {
- // Pivoting
- int maxRow = col;
- for (int row = col + 1; row < n; row++)
- if (Math.Abs(A[row, col]) > Math.Abs(A[maxRow, col]))
- maxRow = row;
- if (maxRow != col)
- {
- for (int j = 0; j < n; j++)
- (A[col, j], A[maxRow, j]) = (A[maxRow, j], A[col, j]);
- (b[col], b[maxRow]) = (b[maxRow], b[col]);
- }
-
- for (int row = col + 1; row < n; row++)
- {
- double factor = A[row, col] / A[col, col];
- for (int j = col; j < n; j++)
- A[row, j] -= factor * A[col, j];
- b[row] -= factor * b[col];
- }
- }
+ // The not-a-knot rows break the tridiagonal structure, so this system is solved
+ // densely. The tridiagonal branch above keeps the Thomas algorithm, which is O(n)
+ // against LU's O(n^3) and must not be replaced by it.
+ var solution = new LuDecomposition(new Matrix(A)).Solve(new VectorN(b));
- // Back substitution
- for (int i = n - 1; i >= 0; i--)
- {
- double sum = b[i];
- for (int j = i + 1; j < n; j++)
- sum -= A[i, j] * M[j];
- M[i] = sum / A[i, i];
- }
+ for (int i = 0; i < n; i++)
+ M[i] = solution[i];
}
private static void SolveTridiagonal(double[] lower, double[] diag, double[] upper, double[] rhs, double[] result, int n)
diff --git a/Numerics/Numerics/Numerics/Interpolation/MultivariateInterpolation.cs b/Numerics/Numerics/Numerics/Interpolation/MultivariateInterpolation.cs
index 200804c..bb68106 100644
--- a/Numerics/Numerics/Numerics/Interpolation/MultivariateInterpolation.cs
+++ b/Numerics/Numerics/Numerics/Interpolation/MultivariateInterpolation.cs
@@ -1,5 +1,7 @@
-using System;
+using System;
using System.Linq;
+using CSharpNumerics.Numerics.LinearAlgebra.Decompositions;
+using CSharpNumerics.Numerics.Objects;
namespace CSharpNumerics.Numerics.Interpolation;
@@ -27,6 +29,11 @@ namespace CSharpNumerics.Numerics.Interpolation;
///
public class MultivariateInterpolation
{
+ ///
+ /// Pivot magnitude below which the RBF matrix is treated as numerically singular.
+ ///
+ private const double SingularPivotTolerance = 1e-14;
+
private readonly double[][] _points;
private readonly double[] _values;
private readonly int _n;
@@ -261,51 +268,25 @@ private static double KernelFunction(double r, double epsilon, RbfKernel kernel)
};
}
+ ///
+ /// Solves A x = b for the RBF weights via LU decomposition with partial pivoting.
+ ///
+ ///
+ /// RBF interpolation matrices are ill-conditioned by nature — a poorly chosen epsilon
+ /// makes them numerically singular — so the smallest pivot is checked before solving
+ /// and the caller is told to change epsilon or kernel rather than handed a meaningless
+ /// set of weights.
+ ///
private static double[] SolveLinearSystem(double[,] A, double[] b)
{
- int n = b.Length;
- var aug = new double[n, n + 1];
- for (int i = 0; i < n; i++)
- {
- for (int j = 0; j < n; j++)
- aug[i, j] = A[i, j];
- aug[i, n] = b[i];
- }
+ var lu = new LuDecomposition(new Matrix(A));
- // Gaussian elimination with partial pivoting
- for (int col = 0; col < n; col++)
+ if (lu.SmallestPivotMagnitude < SingularPivotTolerance)
{
- int maxRow = col;
- for (int row = col + 1; row < n; row++)
- if (Math.Abs(aug[row, col]) > Math.Abs(aug[maxRow, col]))
- maxRow = row;
-
- if (maxRow != col)
- for (int j = 0; j <= n; j++)
- (aug[col, j], aug[maxRow, j]) = (aug[maxRow, j], aug[col, j]);
-
- double pivot = aug[col, col];
- if (Math.Abs(pivot) < 1e-14)
- throw new InvalidOperationException("Singular or near-singular RBF matrix. Try different epsilon or kernel.");
-
- for (int row = col + 1; row < n; row++)
- {
- double factor = aug[row, col] / pivot;
- for (int j = col; j <= n; j++)
- aug[row, j] -= factor * aug[col, j];
- }
+ throw new InvalidOperationException("Singular or near-singular RBF matrix. Try different epsilon or kernel.");
}
- // Back substitution
- var x = new double[n];
- for (int i = n - 1; i >= 0; i--)
- {
- double sum = aug[i, n];
- for (int j = i + 1; j < n; j++)
- sum -= aug[i, j] * x[j];
- x[i] = sum / aug[i, i];
- }
- return x;
+ return lu.Solve(new VectorN(b)).Values;
}
private static int FindGridIndex(double[] grid, double val)
diff --git a/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/LuDecomposition.cs b/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/LuDecomposition.cs
index 0cce65f..c2e0f6b 100644
--- a/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/LuDecomposition.cs
+++ b/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/LuDecomposition.cs
@@ -79,18 +79,34 @@ public LuDecomposition(Matrix matrix)
///
/// True if the matrix is singular (a zero pivot was encountered) and cannot be solved.
///
- public bool IsSingular
+ public bool IsSingular => SmallestPivotMagnitude == 0.0;
+
+ ///
+ /// Magnitude of the smallest pivot on the diagonal of U.
+ ///
+ ///
+ /// Exactly zero means the matrix is singular. A very small but non-zero value means it is
+ /// near-singular: a solve will succeed but divide by that pivot, amplifying any error in the
+ /// right-hand side accordingly. Callers working with matrices that are ill-conditioned by
+ /// nature — radial basis function interpolation, for instance — can test this and refuse
+ /// rather than return a meaningless answer.
+ ///
+ public double SmallestPivotMagnitude
{
get
{
+ var smallest = double.PositiveInfinity;
+
for (var j = 0; j < n; j++)
{
- if (lu[j, j] == 0.0)
+ var magnitude = Math.Abs(lu[j, j]);
+ if (magnitude < smallest)
{
- return true;
+ smallest = magnitude;
}
}
- return false;
+
+ return smallest;
}
}
diff --git a/Numerics/Numerics/Numerics/LinearAlgebra/MatrixDecompositionExtensions.cs b/Numerics/Numerics/Numerics/LinearAlgebra/MatrixDecompositionExtensions.cs
index 554b032..4bcd052 100644
--- a/Numerics/Numerics/Numerics/LinearAlgebra/MatrixDecompositionExtensions.cs
+++ b/Numerics/Numerics/Numerics/LinearAlgebra/MatrixDecompositionExtensions.cs
@@ -28,4 +28,23 @@ public static class MatrixDecompositionExtensions
/// Symmetric matrices give real, ascending eigenvalues with orthonormal eigenvectors.
///
public static EigenDecomposition Eigen(this Matrix matrix) => new EigenDecomposition(matrix);
+
+ ///
+ /// Solves A·X = B for a matrix expected to be symmetric positive definite, one solve per
+ /// column of .
+ ///
+ ///
+ /// Uses Cholesky, which is about twice as fast as LU and preserves symmetry, and falls back
+ /// to LU if the matrix turns out not to be positive definite. The fallback matters for
+ /// quantities that are positive definite in theory but can drift in floating point —
+ /// a recursively updated covariance, for instance.
+ ///
+ public static Matrix SolveSymmetricPositiveDefinite(this Matrix matrix, Matrix rightHandSides)
+ {
+ var cholesky = new CholeskyDecomposition(matrix);
+
+ return cholesky.IsPositiveDefinite
+ ? cholesky.Solve(rightHandSides)
+ : new LuDecomposition(matrix).Solve(rightHandSides);
+ }
}
diff --git a/Numerics/Numerics/Numerics/NumericExtensions.cs b/Numerics/Numerics/Numerics/NumericExtensions.cs
index 0b7bbc6..c8619e4 100644
--- a/Numerics/Numerics/Numerics/NumericExtensions.cs
+++ b/Numerics/Numerics/Numerics/NumericExtensions.cs
@@ -1,5 +1,6 @@
using System.Collections.Generic;
using System;
+using CSharpNumerics.Numerics.RootFinding;
namespace CSharpNumerics.Numerics;
@@ -153,18 +154,14 @@ public static double Dot(this List firstList, List secondList)
/// Function for which to find a root.
/// Initial guess.
/// An approximation of x such that func(x) is near zero.
- public static double NewtonRaphson(this Func func, double xZero = 1.0)
- {
- var value = xZero;
-
- for (var j = 0; j < 100; j++)
- {
- var y = func(value);
- var yPrime = func.Derivate(value);
-
- value -= y / yPrime;
- }
-
- return value;
- }
+ ///
+ /// Delegates to , which
+ /// stops as soon as it converges instead of always running the full iteration budget, and gives up
+ /// where the derivative vanishes instead of iterating on NaN. The returned value is the method's
+ /// best estimate whether or not it converged — use
+ /// directly, or for a bracketed search, when you need
+ /// to know which.
+ ///
+ public static double NewtonRaphson(this Func func, double xZero = 1.0) =>
+ RootFinder.Newton(func, xZero).Value;
}
diff --git a/Numerics/Numerics/Numerics/README.md b/Numerics/Numerics/Numerics/README.md
index ff5260b..2d756ac 100644
--- a/Numerics/Numerics/Numerics/README.md
+++ b/Numerics/Numerics/Numerics/README.md
@@ -14,6 +14,9 @@ Func func = x => Math.Pow(x, 2) - 4;
double root = func.NewtonRaphson(); // 2
```
+See [Root Finding](#-root-finding) for Brent, bisection and the secant method, and for
+checking whether a solve actually converged.
+
**Number Theory**
```csharp
@@ -1051,6 +1054,62 @@ var dominant = A.DominantEigenVector();
var solution = matrix.GaussElimination(vector);
```
+---
+
+## 🎯 Root Finding
+
+Solvers for `f(x) = 0`, in `CSharpNumerics.Numerics.RootFinding`. Every method returns a
+`RootResult` that says whether it actually converged — a failed run gives back its best
+estimate, never a silent `NaN`.
+
+```csharp
+using CSharpNumerics.Numerics.RootFinding;
+
+Func f = x => Math.Cos(x) - x;
+
+// Brent's method — the default when you can bracket the root.
+// Always converges, and nearly as fast as Newton.
+var result = f.FindRoot(0, 1);
+if (result.Converged)
+{
+ double root = result.Value; // 0.7390851332151607
+ int iterations = result.Iterations;
+ double residual = result.Residual; // |f(root)|
+}
+
+// Or throw instead of branching:
+double x = f.FindRoot(0, 1).EnsureConverged();
+```
+
+Pick the method to match what you know about the function:
+
+| Method | Needs | Converges | Use when |
+|--------|-------|-----------|----------|
+| `RootFinder.Brent` | A bracket | Always | The default choice |
+| `RootFinder.Bisection` | A bracket | Always, linearly | The function is nasty and speed does not matter |
+| `RootFinder.Secant` | Two starting points | May diverge | No bracket, no derivative |
+| `RootFinder.Newton` | One starting point | May diverge | A good guess, ideally with an analytic derivative |
+
+```csharp
+var bisection = RootFinder.Bisection(f, 0, 1);
+var secant = RootFinder.Secant(f, 0, 1);
+var newton = RootFinder.Newton(f, initialGuess: 0.5);
+
+// An analytic derivative is both faster and more accurate than the
+// finite-difference default.
+var exact = RootFinder.Newton(f, x => -Math.Sin(x) - 1, 0.5);
+
+// Tolerance and iteration budget are configurable on every method.
+var tight = RootFinder.Brent(f, 0, 1, tolerance: 1e-15, maxIterations: 200);
+```
+
+Newton's method stops and reports failure where the derivative vanishes or the iterate stops
+being finite, rather than dividing into an infinity and iterating on `NaN`. Bracketed methods
+throw if the interval does not contain a sign change.
+
+The older `func.NewtonRaphson()` extension still works unchanged and now delegates here, so it
+exits as soon as it converges instead of always running its full iteration budget.
+
---
## ✨ Interpolation
diff --git a/Numerics/Numerics/Numerics/RootFinding/RootFinder.cs b/Numerics/Numerics/Numerics/RootFinding/RootFinder.cs
new file mode 100644
index 0000000..8f0f3cf
--- /dev/null
+++ b/Numerics/Numerics/Numerics/RootFinding/RootFinder.cs
@@ -0,0 +1,460 @@
+using System;
+
+namespace CSharpNumerics.Numerics.RootFinding;
+
+///
+/// Scalar root finders for f(x) = 0.
+///
+///
+///
+/// Every method reports whether it converged rather than silently returning a bad value —
+/// see .
+///
+///
+/// Choosing a method: is the default choice when a bracket is known.
+/// It cannot fail to converge and is nearly as fast as Newton's method.
+/// is slower but the most robust. and need only a
+/// starting point rather than a bracket, but may diverge.
+///
+///
+public static class RootFinder
+{
+ /// Default convergence tolerance.
+ public const double DefaultTolerance = 1e-12;
+
+ /// Default iteration budget.
+ public const int DefaultMaxIterations = 100;
+
+ ///
+ /// Machine epsilon for (2⁻⁵²).
+ /// Note this is not , which is the smallest subnormal value.
+ ///
+ private const double MachineEpsilon = 2.220446049250313e-16;
+
+ ///
+ /// Step used by the default finite-difference derivative in .
+ /// Matches at order 1.
+ ///
+ private const double DerivativeStep = 10.0 * DerivativeExtensions.h;
+
+ ///
+ /// Bisection: repeatedly halves a bracketing interval. Converges for any continuous
+ /// function that changes sign across the bracket, but only linearly — roughly one
+ /// extra correct bit per iteration.
+ ///
+ /// The function whose root is sought.
+ /// Lower end of a bracketing interval.
+ /// Upper end of a bracketing interval.
+ /// Convergence tolerance on the half-width of the interval.
+ /// Maximum number of iterations.
+ /// is null.
+ /// The interval is degenerate, or does not bracket a sign change.
+ public static RootResult Bisection(
+ Func function,
+ double lower,
+ double upper,
+ double tolerance = DefaultTolerance,
+ int maxIterations = DefaultMaxIterations)
+ {
+ ValidateBracket(function, ref lower, ref upper, tolerance, maxIterations, out var fLower, out var fUpper);
+
+ if (fLower == 0.0)
+ {
+ return new RootResult(lower, true, 0, 0.0, nameof(Bisection));
+ }
+
+ if (fUpper == 0.0)
+ {
+ return new RootResult(upper, true, 0, 0.0, nameof(Bisection));
+ }
+
+ var a = lower;
+ var b = upper;
+ var fa = fLower;
+ var midpoint = a;
+ var fMid = fa;
+
+ for (var iteration = 1; iteration <= maxIterations; iteration++)
+ {
+ midpoint = 0.5 * (a + b);
+ fMid = function(midpoint);
+
+ if (fMid == 0.0 || 0.5 * (b - a) <= tolerance)
+ {
+ return new RootResult(midpoint, true, iteration, Math.Abs(fMid), nameof(Bisection));
+ }
+
+ if (Math.Sign(fMid) == Math.Sign(fa))
+ {
+ a = midpoint;
+ fa = fMid;
+ }
+ else
+ {
+ b = midpoint;
+ }
+ }
+
+ return new RootResult(midpoint, false, maxIterations, Math.Abs(fMid), nameof(Bisection));
+ }
+
+ ///
+ /// Secant method: Newton's method with the derivative replaced by a finite difference
+ /// over the two previous iterates. Needs no derivative and no bracket, but may diverge.
+ ///
+ /// The function whose root is sought.
+ /// First starting point.
+ /// Second starting point, distinct from the first.
+ /// Convergence tolerance on the step size and residual.
+ /// Maximum number of iterations.
+ /// is null.
+ /// The two starting points are equal.
+ public static RootResult Secant(
+ Func function,
+ double first,
+ double second,
+ double tolerance = DefaultTolerance,
+ int maxIterations = DefaultMaxIterations)
+ {
+ ValidateCommon(function, tolerance, maxIterations);
+
+ if (first == second)
+ {
+ throw new ArgumentException("The two starting points must be distinct.", nameof(second));
+ }
+
+ var previous = first;
+ var current = second;
+ var fPrevious = function(previous);
+ var fCurrent = function(current);
+
+ for (var iteration = 1; iteration <= maxIterations; iteration++)
+ {
+ if (Math.Abs(fCurrent) <= tolerance)
+ {
+ return new RootResult(current, true, iteration - 1, Math.Abs(fCurrent), nameof(Secant));
+ }
+
+ var denominator = fCurrent - fPrevious;
+
+ // The secant through two equal function values is horizontal: no step is defined.
+ if (denominator == 0.0)
+ {
+ return new RootResult(current, false, iteration - 1, Math.Abs(fCurrent), nameof(Secant));
+ }
+
+ var step = fCurrent * (current - previous) / denominator;
+ var next = current - step;
+
+ if (!IsFinite(next))
+ {
+ return new RootResult(current, false, iteration - 1, Math.Abs(fCurrent), nameof(Secant));
+ }
+
+ previous = current;
+ fPrevious = fCurrent;
+ current = next;
+ fCurrent = function(current);
+
+ if (Math.Abs(step) <= tolerance * (1.0 + Math.Abs(current)))
+ {
+ return new RootResult(current, true, iteration, Math.Abs(fCurrent), nameof(Secant));
+ }
+ }
+
+ return new RootResult(current, false, maxIterations, Math.Abs(fCurrent), nameof(Secant));
+ }
+
+ ///
+ /// Brent's method: combines bisection, the secant method and inverse quadratic
+ /// interpolation. Keeps the root bracketed at all times, so it always converges,
+ /// while achieving superlinear speed on well-behaved functions. The default choice
+ /// when a bracket is available.
+ ///
+ /// The function whose root is sought.
+ /// Lower end of a bracketing interval.
+ /// Upper end of a bracketing interval.
+ /// Convergence tolerance on the root position.
+ /// Maximum number of iterations.
+ /// is null.
+ /// The interval is degenerate, or does not bracket a sign change.
+ public static RootResult Brent(
+ Func function,
+ double lower,
+ double upper,
+ double tolerance = DefaultTolerance,
+ int maxIterations = DefaultMaxIterations)
+ {
+ ValidateBracket(function, ref lower, ref upper, tolerance, maxIterations, out var fa, out var fb);
+
+ if (fa == 0.0)
+ {
+ return new RootResult(lower, true, 0, 0.0, nameof(Brent));
+ }
+
+ if (fb == 0.0)
+ {
+ return new RootResult(upper, true, 0, 0.0, nameof(Brent));
+ }
+
+ var a = lower;
+ var b = upper;
+ var c = a;
+ var fc = fa;
+ var d = b - a;
+ var e = d;
+
+ for (var iteration = 1; iteration <= maxIterations; iteration++)
+ {
+ // Keep c on the opposite side of the root from b.
+ if (Math.Sign(fb) == Math.Sign(fc))
+ {
+ c = a;
+ fc = fa;
+ d = b - a;
+ e = d;
+ }
+
+ // Keep b as the best estimate so far.
+ if (Math.Abs(fc) < Math.Abs(fb))
+ {
+ a = b;
+ b = c;
+ c = a;
+ fa = fb;
+ fb = fc;
+ fc = fa;
+ }
+
+ var tolerance1 = 2.0 * MachineEpsilon * Math.Abs(b) + 0.5 * tolerance;
+ var bisectionStep = 0.5 * (c - b);
+
+ if (Math.Abs(bisectionStep) <= tolerance1 || fb == 0.0)
+ {
+ return new RootResult(b, true, iteration, Math.Abs(fb), nameof(Brent));
+ }
+
+ if (Math.Abs(e) >= tolerance1 && Math.Abs(fa) > Math.Abs(fb))
+ {
+ // Attempt interpolation: linear when only two points are distinct,
+ // inverse quadratic when three are.
+ double numerator;
+ double denominator;
+ var s = fb / fa;
+
+ if (a == c)
+ {
+ numerator = 2.0 * bisectionStep * s;
+ denominator = 1.0 - s;
+ }
+ else
+ {
+ var q = fa / fc;
+ var r = fb / fc;
+ numerator = s * (2.0 * bisectionStep * q * (q - r) - (b - a) * (r - 1.0));
+ denominator = (q - 1.0) * (r - 1.0) * (s - 1.0);
+ }
+
+ if (numerator > 0.0)
+ {
+ denominator = -denominator;
+ }
+
+ numerator = Math.Abs(numerator);
+
+ // Accept the interpolated step only if it stays inside the bracket and
+ // improves on the step before last; otherwise fall back to bisection.
+ var bound = 3.0 * bisectionStep * denominator - Math.Abs(tolerance1 * denominator);
+ var previousStep = Math.Abs(e * denominator);
+
+ if (2.0 * numerator < Math.Min(bound, previousStep))
+ {
+ e = d;
+ d = numerator / denominator;
+ }
+ else
+ {
+ d = bisectionStep;
+ e = d;
+ }
+ }
+ else
+ {
+ d = bisectionStep;
+ e = d;
+ }
+
+ a = b;
+ fa = fb;
+ b += Math.Abs(d) > tolerance1
+ ? d
+ : bisectionStep >= 0.0 ? tolerance1 : -tolerance1;
+ fb = function(b);
+ }
+
+ return new RootResult(b, false, maxIterations, Math.Abs(fb), nameof(Brent));
+ }
+
+ ///
+ /// Newton's method using a finite-difference derivative. Converges quadratically near
+ /// a simple root, but can diverge from a poor starting point and stalls where the
+ /// derivative vanishes — both are reported rather than hidden.
+ ///
+ /// The function whose root is sought.
+ /// Starting point.
+ /// Convergence tolerance on the step size and residual.
+ /// Maximum number of iterations.
+ ///
+ /// The derivative is approximated by the same backward difference that
+ /// uses at
+ /// order 1, so results match. Prefer the overload taking an analytic derivative when one
+ /// is available: it is both faster and more accurate.
+ ///
+ /// is null.
+ public static RootResult Newton(
+ Func function,
+ double initialGuess = 1.0,
+ double tolerance = DefaultTolerance,
+ int maxIterations = DefaultMaxIterations)
+ {
+ if (function == null)
+ {
+ throw new ArgumentNullException(nameof(function));
+ }
+
+ return Newton(
+ function,
+ x => (function(x) - function(x - DerivativeStep)) / DerivativeStep,
+ initialGuess,
+ tolerance,
+ maxIterations);
+ }
+
+ ///
+ /// Newton's method with an analytic derivative.
+ ///
+ /// The function whose root is sought.
+ /// The derivative of .
+ /// Starting point.
+ /// Convergence tolerance on the step size and residual.
+ /// Maximum number of iterations.
+ /// or is null.
+ public static RootResult Newton(
+ Func function,
+ Func derivative,
+ double initialGuess,
+ double tolerance = DefaultTolerance,
+ int maxIterations = DefaultMaxIterations)
+ {
+ ValidateCommon(function, tolerance, maxIterations);
+
+ if (derivative == null)
+ {
+ throw new ArgumentNullException(nameof(derivative));
+ }
+
+ var x = initialGuess;
+ var fx = function(x);
+
+ for (var iteration = 1; iteration <= maxIterations; iteration++)
+ {
+ if (Math.Abs(fx) <= tolerance)
+ {
+ return new RootResult(x, true, iteration - 1, Math.Abs(fx), nameof(Newton));
+ }
+
+ var slope = derivative(x);
+
+ // A vanishing or non-finite slope gives no usable step — stop instead of
+ // dividing into an infinity and iterating on NaN.
+ if (slope == 0.0 || !IsFinite(slope))
+ {
+ return new RootResult(x, false, iteration - 1, Math.Abs(fx), nameof(Newton));
+ }
+
+ var step = fx / slope;
+ var next = x - step;
+
+ if (!IsFinite(next))
+ {
+ return new RootResult(x, false, iteration - 1, Math.Abs(fx), nameof(Newton));
+ }
+
+ x = next;
+ fx = function(x);
+
+ if (Math.Abs(step) <= tolerance * (1.0 + Math.Abs(x)))
+ {
+ return new RootResult(x, true, iteration, Math.Abs(fx), nameof(Newton));
+ }
+ }
+
+ return new RootResult(x, false, maxIterations, Math.Abs(fx), nameof(Newton));
+ }
+
+ private static void ValidateCommon(Func function, double tolerance, int maxIterations)
+ {
+ if (function == null)
+ {
+ throw new ArgumentNullException(nameof(function));
+ }
+
+ if (tolerance <= 0.0 || !IsFinite(tolerance))
+ {
+ throw new ArgumentOutOfRangeException(nameof(tolerance), "Tolerance must be positive and finite.");
+ }
+
+ if (maxIterations < 1)
+ {
+ throw new ArgumentOutOfRangeException(nameof(maxIterations), "At least one iteration is required.");
+ }
+ }
+
+ private static void ValidateBracket(
+ Func function,
+ ref double lower,
+ ref double upper,
+ double tolerance,
+ int maxIterations,
+ out double fLower,
+ out double fUpper)
+ {
+ ValidateCommon(function, tolerance, maxIterations);
+
+ if (!IsFinite(lower) || !IsFinite(upper))
+ {
+ throw new ArgumentException("The bracket endpoints must be finite.");
+ }
+
+ if (lower == upper)
+ {
+ throw new ArgumentException("The bracket endpoints must be distinct.");
+ }
+
+ if (lower > upper)
+ {
+ var swap = lower;
+ lower = upper;
+ upper = swap;
+ }
+
+ fLower = function(lower);
+ fUpper = function(upper);
+
+ if (!IsFinite(fLower) || !IsFinite(fUpper))
+ {
+ throw new ArgumentException("The function must be finite at both bracket endpoints.");
+ }
+
+ // A sign change guarantees a root for a continuous function; without one these
+ // methods have nothing to bisect towards.
+ if (fLower != 0.0 && fUpper != 0.0 && Math.Sign(fLower) == Math.Sign(fUpper))
+ {
+ throw new ArgumentException(
+ $"The interval [{lower:G6}, {upper:G6}] does not bracket a sign change: " +
+ $"f(lower) = {fLower:G6} and f(upper) = {fUpper:G6} have the same sign.");
+ }
+ }
+
+ private static bool IsFinite(double value) => !double.IsNaN(value) && !double.IsInfinity(value);
+}
diff --git a/Numerics/Numerics/Numerics/RootFinding/RootFindingExtensions.cs b/Numerics/Numerics/Numerics/RootFinding/RootFindingExtensions.cs
new file mode 100644
index 0000000..d3873c8
--- /dev/null
+++ b/Numerics/Numerics/Numerics/RootFinding/RootFindingExtensions.cs
@@ -0,0 +1,34 @@
+using System;
+
+namespace CSharpNumerics.Numerics.RootFinding;
+
+///
+/// Extension-method facade over , in the same style as the other
+/// numerics extensions.
+///
+public static class RootFindingExtensions
+{
+ ///
+ /// Finds a root of the function in the bracketing interval [,
+ /// ] using Brent's method.
+ ///
+ /// The function whose root is sought.
+ /// Lower end of a bracketing interval.
+ /// Upper end of a bracketing interval.
+ /// Convergence tolerance on the root position.
+ /// Maximum number of iterations.
+ /// The outcome, including whether the run converged.
+ ///
+ ///
+ /// Func<double, double> f = x => x * x - 4;
+ /// var root = f.FindRoot(0, 5); // root.Value ≈ 2
+ ///
+ ///
+ public static RootResult FindRoot(
+ this Func function,
+ double lower,
+ double upper,
+ double tolerance = RootFinder.DefaultTolerance,
+ int maxIterations = RootFinder.DefaultMaxIterations) =>
+ RootFinder.Brent(function, lower, upper, tolerance, maxIterations);
+}
diff --git a/Numerics/Numerics/Numerics/RootFinding/RootResult.cs b/Numerics/Numerics/Numerics/RootFinding/RootResult.cs
new file mode 100644
index 0000000..fa86e2f
--- /dev/null
+++ b/Numerics/Numerics/Numerics/RootFinding/RootResult.cs
@@ -0,0 +1,61 @@
+using System;
+
+namespace CSharpNumerics.Numerics.RootFinding;
+
+///
+/// Outcome of a root-finding run: the approximate root together with the information
+/// needed to judge whether it can be trusted.
+///
+///
+/// A root finder that stops without converging still returns its best estimate, so the
+/// value alone says nothing about success — always check before
+/// using it. Use for the raw estimate or
+/// to turn a failure into an exception.
+///
+public sealed class RootResult
+{
+ public RootResult(double value, bool converged, int iterations, double residual, string method)
+ {
+ Value = value;
+ Converged = converged;
+ Iterations = iterations;
+ Residual = residual;
+ Method = method;
+ }
+
+ /// The approximate root. Only meaningful when is true.
+ public double Value { get; }
+
+ /// True if the method reached its tolerance before exhausting its iteration budget.
+ public bool Converged { get; }
+
+ /// Number of iterations actually performed.
+ public int Iterations { get; }
+
+ /// The residual |f()| at the returned point.
+ public double Residual { get; }
+
+ /// Name of the method that produced the result, e.g. "Brent".
+ public string Method { get; }
+
+ ///
+ /// Returns if the run converged, otherwise throws.
+ /// Use when a non-converged result is a programming error rather than something to handle.
+ ///
+ /// The run did not converge.
+ public double EnsureConverged()
+ {
+ if (!Converged)
+ {
+ throw new InvalidOperationException(
+ $"{Method} did not converge after {Iterations} iterations (residual {Residual:G6}).");
+ }
+
+ return Value;
+ }
+
+ public override string ToString() =>
+ Converged
+ ? $"{Method}: x = {Value:G10} after {Iterations} iterations (residual {Residual:G6})"
+ : $"{Method}: did not converge after {Iterations} iterations (best x = {Value:G10}, residual {Residual:G6})";
+}
diff --git a/Numerics/Numerics/Physics/Astro/KeplerOrbit.cs b/Numerics/Numerics/Physics/Astro/KeplerOrbit.cs
index 89ceac9..4c4aa43 100644
--- a/Numerics/Numerics/Physics/Astro/KeplerOrbit.cs
+++ b/Numerics/Numerics/Physics/Astro/KeplerOrbit.cs
@@ -1,4 +1,5 @@
using System;
+using CSharpNumerics.Numerics.RootFinding;
using CSharpNumerics.Physics.Constants;
namespace CSharpNumerics.Physics.Astro;
@@ -19,15 +20,14 @@ public static double TrueAnomaly(double meanAnomaly, double eccentricity, double
if (eccentricity < 0 || eccentricity >= 1)
throw new ArgumentOutOfRangeException(nameof(eccentricity), "Eccentricity must be in [0, 1).");
- // Solve Kepler's equation: M = E - e·sin(E)
- double E = meanAnomaly; // initial guess
- for (int i = 0; i < maxIterations; i++)
- {
- double dE = (E - eccentricity * Math.Sin(E) - meanAnomaly) / (1.0 - eccentricity * Math.Cos(E));
- E -= dE;
- if (Math.Abs(dE) < tolerance)
- break;
- }
+ // Solve Kepler's equation M = E - e·sin(E) by Newton's method. The derivative
+ // 1 - e·cos(E) is bounded below by 1 - e > 0, so no step can stall.
+ double E = RootFinder.Newton(
+ eccentric => eccentric - eccentricity * Math.Sin(eccentric) - meanAnomaly,
+ eccentric => 1.0 - eccentricity * Math.Cos(eccentric),
+ meanAnomaly,
+ tolerance,
+ maxIterations).Value;
// Convert eccentric anomaly to true anomaly
double sinNu = Math.Sqrt(1.0 - eccentricity * eccentricity) * Math.Sin(E) / (1.0 - eccentricity * Math.Cos(E));
diff --git a/Numerics/Numerics/Physics/Gravitation/LagrangePoints.cs b/Numerics/Numerics/Physics/Gravitation/LagrangePoints.cs
index c31453a..b0cd9f2 100644
--- a/Numerics/Numerics/Physics/Gravitation/LagrangePoints.cs
+++ b/Numerics/Numerics/Physics/Gravitation/LagrangePoints.cs
@@ -1,5 +1,6 @@
using System;
using CSharpNumerics.Numerics.Objects;
+using CSharpNumerics.Numerics.RootFinding;
namespace CSharpNumerics.Physics.Gravitation;
@@ -58,31 +59,30 @@ public static (Vector L1, Vector L2, Vector L3, Vector L4, Vector L5) All(
// Newton's method on f(x) = x − (1−mu)(x+mu)/|x+mu|³ − mu(x−1+mu)/|x−1+mu|³.
private static double SolveCollinear(double initialGuess, double mu)
{
- double x = initialGuess;
double m1Pos = -mu; // larger mass position
double m2Pos = 1.0 - mu; // smaller mass position
- for (int iter = 0; iter < 100; iter++)
+ double Equation(double x)
{
- double d1 = x - m1Pos;
- double d2 = x - m2Pos;
- double a1 = Math.Abs(d1);
- double a2 = Math.Abs(d2);
+ double a1 = Math.Abs(x - m1Pos);
+ double a2 = Math.Abs(x - m2Pos);
- double f = x
- - (1.0 - mu) * d1 / (a1 * a1 * a1)
- - mu * d2 / (a2 * a2 * a2);
+ return x
+ - (1.0 - mu) * (x - m1Pos) / (a1 * a1 * a1)
+ - mu * (x - m2Pos) / (a2 * a2 * a2);
+ }
- // f'(x) = 1 + 2(1−mu)/|d1|³ + 2·mu/|d2|³ (always positive)
- double fp = 1.0
- + 2.0 * (1.0 - mu) / (a1 * a1 * a1)
- + 2.0 * mu / (a2 * a2 * a2);
+ // f'(x) = 1 + 2(1−mu)/|d1|³ + 2·mu/|d2|³ (always positive)
+ double Derivative(double x)
+ {
+ double a1 = Math.Abs(x - m1Pos);
+ double a2 = Math.Abs(x - m2Pos);
- double dx = f / fp;
- x -= dx;
- if (Math.Abs(dx) < 1e-14) break;
+ return 1.0
+ + 2.0 * (1.0 - mu) / (a1 * a1 * a1)
+ + 2.0 * mu / (a2 * a2 * a2);
}
- return x;
+ return RootFinder.Newton(Equation, Derivative, initialGuess, 1e-14, 100).Value;
}
}
diff --git a/Numerics/Numerics/Physics/Mechanics/Oscillations/CoupledOscillators.cs b/Numerics/Numerics/Physics/Mechanics/Oscillations/CoupledOscillators.cs
index 8bc94f2..9c7eb65 100644
--- a/Numerics/Numerics/Physics/Mechanics/Oscillations/CoupledOscillators.cs
+++ b/Numerics/Numerics/Physics/Mechanics/Oscillations/CoupledOscillators.cs
@@ -1,4 +1,5 @@
-using CSharpNumerics.Numerics;
+using CSharpNumerics.Numerics;
+using CSharpNumerics.Numerics.LinearAlgebra.Decompositions;
using CSharpNumerics.Numerics.Objects;
using CSharpNumerics.Statistics.Data;
using System;
@@ -16,7 +17,7 @@ namespace CSharpNumerics.Physics.Mechanics.Oscillations
/// damping matrix.
///
///
- /// Normal modes are computed via the Jacobi eigenvalue algorithm on the
+ /// Normal modes are computed with EigenDecomposition on the
/// symmetrised dynamical matrix D = L⁻¹KL⁻¹ where L = diag(√m_i),
/// giving eigenvalues ω² and orthonormal mode shapes.
///
@@ -248,7 +249,7 @@ public Matrix MassMatrix()
#endregion
- #region Eigen-decomposition (Jacobi)
+ #region Eigen-decomposition
///
/// Builds the symmetrised dynamical matrix D = L⁻¹ K L⁻¹
@@ -265,110 +266,6 @@ public Matrix MassMatrix()
return D;
}
- ///
- /// Jacobi eigenvalue algorithm for real symmetric matrices.
- /// Returns eigenvalues in ascending order and orthonormal eigenvectors
- /// as columns of the eigenvectors matrix.
- ///
- private void JacobiEigen(double[,] A, int n, out double[] eigenvalues, out double[,] eigenvectors)
- {
- // Work on a copy
- var S = new double[n, n];
- Array.Copy(A, S, A.Length);
-
- // Eigenvector accumulator (starts as identity)
- var V = new double[n, n];
- for (int i = 0; i < n; i++)
- V[i, i] = 1.0;
-
- int maxIterations = 100 * n * n;
- double tol = 1e-12;
-
- for (int iter = 0; iter < maxIterations; iter++)
- {
- // Find largest off-diagonal element
- int p = 0, q = 1;
- double maxVal = 0;
- for (int i = 0; i < n; i++)
- {
- for (int j = i + 1; j < n; j++)
- {
- double absVal = Math.Abs(S[i, j]);
- if (absVal > maxVal)
- {
- maxVal = absVal;
- p = i;
- q = j;
- }
- }
- }
-
- if (maxVal < tol) break;
-
- // Compute rotation angle
- double theta;
- if (Math.Abs(S[p, p] - S[q, q]) < 1e-15)
- {
- theta = Math.PI / 4.0;
- }
- else
- {
- theta = 0.5 * Math.Atan2(2.0 * S[p, q], S[p, p] - S[q, q]);
- }
-
- double c = Math.Cos(theta);
- double s = Math.Sin(theta);
-
- // Apply Jacobi rotation: S' = G^T S G
- // Update rows/columns p and q
- var Sp = new double[n];
- var Sq = new double[n];
- for (int i = 0; i < n; i++)
- {
- Sp[i] = c * S[p, i] + s * S[q, i];
- Sq[i] = -s * S[p, i] + c * S[q, i];
- }
- for (int i = 0; i < n; i++)
- {
- S[p, i] = Sp[i];
- S[q, i] = Sq[i];
- S[i, p] = Sp[i];
- S[i, q] = Sq[i];
- }
- // Fix the 2x2 block
- double Spp = c * Sp[p] + s * Sp[q];
- double Sqq = -s * Sq[p] + c * Sq[q];
- double Spq = -s * Sp[p] + c * Sp[q];
- S[p, p] = Spp;
- S[q, q] = Sqq;
- S[p, q] = 0;
- S[q, p] = 0;
-
- // Accumulate eigenvectors: V = V * G
- for (int i = 0; i < n; i++)
- {
- double vip = V[i, p];
- double viq = V[i, q];
- V[i, p] = c * vip + s * viq;
- V[i, q] = -s * vip + c * viq;
- }
- }
-
- // Extract eigenvalues and sort by ascending order
- var indices = Enumerable.Range(0, n)
- .OrderBy(i => S[i, i])
- .ToArray();
-
- eigenvalues = new double[n];
- eigenvectors = new double[n, n];
- for (int k = 0; k < n; k++)
- {
- eigenvalues[k] = S[indices[k], indices[k]];
- for (int i = 0; i < n; i++)
- eigenvectors[i, k] = V[i, indices[k]];
- }
- }
-
///
/// Ensures the eigendecomposition is computed and cached.
///
@@ -377,7 +274,13 @@ private void EnsureEigendecomposition()
if (_cachedEigenvalues != null) return;
var D = SymmetricDynamicalMatrix();
- JacobiEigen(D, _n, out var rawEigenvalues, out var rawEigenvectors);
+
+ // D is symmetric by construction, so EigenDecomposition takes its symmetric path:
+ // real eigenvalues in ascending order with orthonormal eigenvectors as columns —
+ // the same contract the private Jacobi solver provided.
+ var eigen = new EigenDecomposition(new Matrix(D));
+ var rawEigenvalues = eigen.RealEigenvalues;
+ var rawEigenvectors = eigen.EigenVectors.values;
// rawEigenvectors are for the symmetrised problem D.
// Physical mode shapes: φ_i = v_i / √m_i
diff --git a/Numerics/Numerics/Statistics/Fitting/FittingSolver.cs b/Numerics/Numerics/Statistics/Fitting/FittingSolver.cs
index 5c6ae32..92c5e82 100644
--- a/Numerics/Numerics/Statistics/Fitting/FittingSolver.cs
+++ b/Numerics/Numerics/Statistics/Fitting/FittingSolver.cs
@@ -1,3 +1,4 @@
+using CSharpNumerics.Numerics.LinearAlgebra.Decompositions;
using CSharpNumerics.Numerics.Objects;
using System;
@@ -68,6 +69,110 @@ internal static double[] Solve(double[,] A, double[] b)
return x;
}
+ ///
+ /// Solves the least squares problem min ‖Xβ − y‖ by QR decomposition, returning the
+ /// coefficients together with (XᵀX)⁻¹ for standard errors.
+ ///
+ ///
+ /// Works on the design matrix directly rather than forming the normal equations XᵀX, which
+ /// square the condition number of X. For a Vandermonde design — any polynomial fit — that
+ /// squaring is a real loss of accuracy. (XᵀX)⁻¹ is recovered as R⁻¹R⁻ᵀ, since XᵀX = RᵀR,
+ /// so the Gram matrix is never formed or inverted at all.
+ ///
+ /// Design matrix (n × p), n > p.
+ /// Response vector (length n).
+ /// Number of parameters.
+ /// The design matrix is rank deficient.
+ internal static (double[] Beta, double[,] GramInverse) SolveLeastSquares(double[,] X, double[] y, int p)
+ {
+ var qr = new QrDecomposition(new Matrix(X));
+
+ if (!qr.IsFullRank)
+ {
+ throw new InvalidOperationException("Singular matrix encountered during fitting.");
+ }
+
+ var beta = qr.Solve(new VectorN(y)).Values;
+
+ return (beta, GramInverseFromR(qr.R, p));
+ }
+
+ ///
+ /// Solves the weighted least squares problem min ‖√W(Xβ − y)‖ by QR, returning the
+ /// coefficients together with (XᵀWX)⁻¹.
+ ///
+ ///
+ /// Scaling each row i by √wᵢ turns the weighted problem into an ordinary one, so the same
+ /// QR path applies and the weighted normal equations are likewise never formed.
+ ///
+ internal static (double[] Beta, double[,] GramInverse) SolveWeightedLeastSquares(
+ double[,] X, double[] y, double[] weights, int n, int p)
+ {
+ var scaledX = new double[n, p];
+ var scaledY = new double[n];
+
+ for (int i = 0; i < n; i++)
+ {
+ // A negative weight has no meaning here and would make the square root undefined.
+ double rootWeight = Math.Sqrt(weights[i] > 0 ? weights[i] : 0.0);
+
+ for (int j = 0; j < p; j++)
+ scaledX[i, j] = rootWeight * X[i, j];
+
+ scaledY[i] = rootWeight * y[i];
+ }
+
+ return SolveLeastSquares(scaledX, scaledY, p);
+ }
+
+ ///
+ /// Computes (XᵀX)⁻¹ = R⁻¹R⁻ᵀ from the triangular factor of X, without forming XᵀX.
+ ///
+ internal static double[,] GramInverseFromR(Matrix r, int p)
+ {
+ // Invert the upper triangular R by back substitution, one column at a time.
+ var rInverse = new double[p, p];
+
+ for (int col = p - 1; col >= 0; col--)
+ {
+ double diagonal = r.values[col, col];
+
+ if (Math.Abs(diagonal) < 1e-14)
+ {
+ throw new InvalidOperationException("Singular matrix cannot be inverted.");
+ }
+
+ rInverse[col, col] = 1.0 / diagonal;
+
+ for (int row = col - 1; row >= 0; row--)
+ {
+ double sum = 0;
+ for (int k = row + 1; k <= col; k++)
+ sum += r.values[row, k] * rInverse[k, col];
+
+ rInverse[row, col] = -sum / r.values[row, row];
+ }
+ }
+
+ // (XᵀX)⁻¹ = R⁻¹ · R⁻ᵀ, symmetric by construction.
+ var gramInverse = new double[p, p];
+
+ for (int i = 0; i < p; i++)
+ {
+ for (int j = i; j < p; j++)
+ {
+ double sum = 0;
+ for (int k = 0; k < p; k++)
+ sum += rInverse[i, k] * rInverse[j, k];
+
+ gramInverse[i, j] = sum;
+ gramInverse[j, i] = sum;
+ }
+ }
+
+ return gramInverse;
+ }
+
/// Inverts a matrix using Gauss-Jordan elimination with partial pivoting.
internal static double[,] Invert(double[,] A, int n)
{
diff --git a/Numerics/Numerics/Statistics/Fitting/LeastSquaresFitter.cs b/Numerics/Numerics/Statistics/Fitting/LeastSquaresFitter.cs
index c522e30..5be7d8e 100644
--- a/Numerics/Numerics/Statistics/Fitting/LeastSquaresFitter.cs
+++ b/Numerics/Numerics/Statistics/Fitting/LeastSquaresFitter.cs
@@ -1,4 +1,4 @@
-using CSharpNumerics.Numerics.Objects;
+using CSharpNumerics.Numerics.Objects;
using System;
namespace CSharpNumerics.Statistics.Fitting;
@@ -47,9 +47,8 @@ public static FittingResult Fit(double[,] designMatrix, VectorN y)
private static FittingResult FitDesignMatrix(double[,] X, double[] y, int n, int p)
{
- double[,] XtX = FittingSolver.MultiplyATA(X, n, p);
- double[] Xty = FittingSolver.MultiplyATb(X, y, n, p);
- double[] beta = FittingSolver.Solve(XtX, Xty);
+ // QR on the design matrix, rather than forming and solving the normal equations.
+ var (beta, XtXInv) = FittingSolver.SolveLeastSquares(X, y, p);
double[] fitted = FittingSolver.ComputeFitted(X, beta, n, p);
double ssRes = 0;
@@ -60,7 +59,6 @@ private static FittingResult FitDesignMatrix(double[,] X, double[] y, int n, int
}
double s2 = (n > p) ? ssRes / (n - p) : 0.0;
- double[,] XtXInv = FittingSolver.Invert(XtX, p);
double[] se = FittingSolver.ComputeStandardErrors(XtXInv, s2, p);
return new FittingResult(
diff --git a/Numerics/Numerics/Statistics/Fitting/NonlinearLeastSquaresFitter.cs b/Numerics/Numerics/Statistics/Fitting/NonlinearLeastSquaresFitter.cs
index dd7e8d6..9d040ba 100644
--- a/Numerics/Numerics/Statistics/Fitting/NonlinearLeastSquaresFitter.cs
+++ b/Numerics/Numerics/Statistics/Fitting/NonlinearLeastSquaresFitter.cs
@@ -1,3 +1,4 @@
+using CSharpNumerics.Numerics.LinearAlgebra.Decompositions;
using CSharpNumerics.Numerics.Objects;
using System;
@@ -165,8 +166,12 @@ private static double[] ComputeStandardErrorsAtSolution(
}
double s2 = (n > p) ? ssRes / (n - p) : 0.0;
- double[,] JtJ = FittingSolver.MultiplyATA(J, n, p);
- double[,] cov = FittingSolver.Invert(JtJ, p);
+ // (JᵀJ)⁻¹ from the triangular factor of J, without forming JᵀJ. The
+ // Levenberg-Marquardt iteration above still works on the damped normal equations:
+ // Marquardt damping is applied to JᵀJ's diagonal, which has no design-matrix
+ // equivalent short of the augmented formulation [J; √(λD)] — a change to the
+ // algorithm rather than to the solver.
+ double[,] cov = FittingSolver.GramInverseFromR(new QrDecomposition(new Matrix(J)).R, p);
return FittingSolver.ComputeStandardErrors(cov, s2, p);
}
catch (InvalidOperationException)
diff --git a/Numerics/Numerics/Statistics/Fitting/ParameterEstimation.cs b/Numerics/Numerics/Statistics/Fitting/ParameterEstimation.cs
index b24ad74..5cb5178 100644
--- a/Numerics/Numerics/Statistics/Fitting/ParameterEstimation.cs
+++ b/Numerics/Numerics/Statistics/Fitting/ParameterEstimation.cs
@@ -1,3 +1,4 @@
+using CSharpNumerics.Numerics.LinearAlgebra.Decompositions;
using CSharpNumerics.Numerics.Objects;
using System;
@@ -56,8 +57,9 @@ public static (double Lower, double Upper) PredictionInterval(
int n = designMatrix.GetLength(0);
int p = designMatrix.GetLength(1);
- double[,] XtX = FittingSolver.MultiplyATA(designMatrix, n, p);
- double[,] XtXInv = FittingSolver.Invert(XtX, p);
+ // (XᵀX)⁻¹ from the triangular factor of X, without forming XᵀX.
+ double[,] XtXInv = FittingSolver.GramInverseFromR(
+ new QrDecomposition(new Matrix(designMatrix)).R, p);
// Residual variance (unbiased)
double ssRes = 0;
diff --git a/Numerics/Numerics/Statistics/Fitting/RobustFitter.cs b/Numerics/Numerics/Statistics/Fitting/RobustFitter.cs
index 81bc755..12a3775 100644
--- a/Numerics/Numerics/Statistics/Fitting/RobustFitter.cs
+++ b/Numerics/Numerics/Statistics/Fitting/RobustFitter.cs
@@ -1,4 +1,4 @@
-using CSharpNumerics.Numerics.Objects;
+using CSharpNumerics.Numerics.Objects;
using System;
namespace CSharpNumerics.Statistics.Fitting;
@@ -76,10 +76,8 @@ private static FittingResult FitDesignMatrix(
if (c <= 0)
c = DefaultC(wf);
- // Step 1: initial OLS fit
- double[,] XtX = FittingSolver.MultiplyATA(X, n, p);
- double[] Xty = FittingSolver.MultiplyATb(X, y, n, p);
- double[] beta = FittingSolver.Solve(XtX, Xty);
+ // Step 1: initial OLS fit, via QR on the design matrix.
+ double[] beta = FittingSolver.SolveLeastSquares(X, y, p).Beta;
double[] weights = new double[n];
@@ -105,12 +103,9 @@ private static FittingResult FitDesignMatrix(
weights[i] = ComputeWeight(u, c, wf);
}
- // WLS step
- double[,] XtWX = FittingSolver.MultiplyATWA(X, weights, n, p);
- double[] XtWy = FittingSolver.MultiplyATWb(X, weights, y, n, p);
-
+ // WLS step, via QR on the row-scaled design matrix.
double[] betaNew;
- try { betaNew = FittingSolver.Solve(XtWX, XtWy); }
+ try { betaNew = FittingSolver.SolveWeightedLeastSquares(X, y, weights, n, p).Beta; }
catch (InvalidOperationException) { break; }
// Check convergence
@@ -136,8 +131,8 @@ private static FittingResult FitDesignMatrix(
double[] se;
try
{
- double[,] XtWX = FittingSolver.MultiplyATWA(X, weights, n, p);
- double[,] XtWXInv = FittingSolver.Invert(XtWX, p);
+ double[,] XtWXInv = FittingSolver
+ .SolveWeightedLeastSquares(X, y, weights, n, p).GramInverse;
se = FittingSolver.ComputeStandardErrors(XtWXInv, s2, p);
}
catch (InvalidOperationException)
diff --git a/Numerics/Numerics/Statistics/Fitting/WeightedLeastSquaresFitter.cs b/Numerics/Numerics/Statistics/Fitting/WeightedLeastSquaresFitter.cs
index 4d11cdb..82ed0e0 100644
--- a/Numerics/Numerics/Statistics/Fitting/WeightedLeastSquaresFitter.cs
+++ b/Numerics/Numerics/Statistics/Fitting/WeightedLeastSquaresFitter.cs
@@ -1,4 +1,4 @@
-using CSharpNumerics.Numerics.Objects;
+using CSharpNumerics.Numerics.Objects;
using System;
namespace CSharpNumerics.Statistics.Fitting;
@@ -44,9 +44,8 @@ public static FittingResult Fit(double[,] designMatrix, VectorN y, VectorN weigh
internal static FittingResult FitDesignMatrix(double[,] X, double[] y, double[] weights, int n, int p)
{
- double[,] XtWX = FittingSolver.MultiplyATWA(X, weights, n, p);
- double[] XtWy = FittingSolver.MultiplyATWb(X, weights, y, n, p);
- double[] beta = FittingSolver.Solve(XtWX, XtWy);
+ // QR on the row-scaled design matrix, rather than the weighted normal equations.
+ var (beta, XtWXInv) = FittingSolver.SolveWeightedLeastSquares(X, y, weights, n, p);
double[] fitted = FittingSolver.ComputeFitted(X, beta, n, p);
double ssRes = 0;
@@ -57,7 +56,6 @@ internal static FittingResult FitDesignMatrix(double[,] X, double[] y, double[]
}
double s2 = (n > p) ? ssRes / (n - p) : 0.0;
- double[,] XtWXInv = FittingSolver.Invert(XtWX, p);
double[] se = FittingSolver.ComputeStandardErrors(XtWXInv, s2, p);
return new FittingResult(
diff --git a/Numerics/Numerics/Statistics/InferentialStatisticsExtensions.cs b/Numerics/Numerics/Statistics/InferentialStatisticsExtensions.cs
index 645bc96..70f4580 100644
--- a/Numerics/Numerics/Statistics/InferentialStatisticsExtensions.cs
+++ b/Numerics/Numerics/Statistics/InferentialStatisticsExtensions.cs
@@ -1,4 +1,6 @@
-using CSharpNumerics.Statistics.Distributions;
+using CSharpNumerics.Numerics.LinearAlgebra.Decompositions;
+using CSharpNumerics.Numerics.Objects;
+using CSharpNumerics.Statistics.Distributions;
using System;
using System.Collections.Generic;
using System.Linq;
@@ -401,46 +403,30 @@ public static StatisticalTestResult OneWayAnova(this IEnumerable enumerabl
// ──────────────────────────────────────────────
///
- /// Solves an augmented matrix [A|b] via Gaussian elimination with partial pivoting.
+ /// Solves an augmented matrix [A|b] via LU decomposition with partial pivoting.
///
+ ///
+ /// Note that the caller forms the normal equations XᵀX before getting here, which squares
+ /// the condition number of the design matrix. Switching the solve to LU removes a duplicate
+ /// elimination routine but does not recover that lost accuracy — only building the design
+ /// matrix and solving it with QR would.
+ ///
private static double[] GaussianElimination(double[,] matrix, int n)
{
- for (int col = 0; col < n; col++)
- {
- // Partial pivot
- int maxRow = col;
- for (int row = col + 1; row < n; row++)
- {
- if (Math.Abs(matrix[row, col]) > Math.Abs(matrix[maxRow, col]))
- maxRow = row;
- }
- for (int j = col; j <= n; j++)
- {
- double tmp = matrix[col, j];
- matrix[col, j] = matrix[maxRow, j];
- matrix[maxRow, j] = tmp;
- }
+ var a = new double[n, n];
+ var b = new double[n];
- // Eliminate below
- for (int row = col + 1; row < n; row++)
+ for (var i = 0; i < n; i++)
+ {
+ for (var j = 0; j < n; j++)
{
- double factor = matrix[row, col] / matrix[col, col];
- for (int j = col; j <= n; j++)
- matrix[row, j] -= factor * matrix[col, j];
+ a[i, j] = matrix[i, j];
}
- }
- // Back substitute
- double[] result = new double[n];
- for (int i = n - 1; i >= 0; i--)
- {
- result[i] = matrix[i, n];
- for (int j = i + 1; j < n; j++)
- result[i] -= matrix[i, j] * result[j];
- result[i] /= matrix[i, i];
+ b[i] = matrix[i, n];
}
- return result;
+ return new LuDecomposition(new Matrix(a)).Solve(new VectorN(b)).Values;
}
///
diff --git a/Numerics/Numerics/Statistics/README.md b/Numerics/Numerics/Statistics/README.md
index 8066b7b..11c32f7 100644
--- a/Numerics/Numerics/Statistics/README.md
+++ b/Numerics/Numerics/Statistics/README.md
@@ -241,7 +241,17 @@ The `CSharpNumerics.Statistics.Fitting` namespace provides a comprehensive curve
### Least Squares (OLS)
-Polynomial and multiple-regression fitting via the normal equations.
+Polynomial and multiple-regression fitting, solved by QR decomposition of the design matrix.
+
+> Ordinary, weighted and robust fitting all solve the design matrix directly rather than forming
+> the normal equations XᵀX, which square its condition number. On a degree-5 polynomial fit over
+> [1, 2] that is worth close to six decimal digits: the normal equations recover the coefficients
+> to 5.6e-5, QR to 7.2e-11. Standard errors come from the triangular factor as R⁻¹R⁻ᵀ, so the Gram
+> matrix is never formed or inverted.
+>
+> The Levenberg–Marquardt iteration in `NonlinearLeastSquaresFitter` is the exception: its damping
+> is applied to JᵀJ's diagonal, so it still works on the damped normal equations. Only its
+> covariance estimate uses QR.
```csharp
using CSharpNumerics.Statistics.Fitting;
diff --git a/Numerics/Numerics/Statistics/StateEstimation/ExtendedKalmanFilter.cs b/Numerics/Numerics/Statistics/StateEstimation/ExtendedKalmanFilter.cs
index 6bebecc..de81b21 100644
--- a/Numerics/Numerics/Statistics/StateEstimation/ExtendedKalmanFilter.cs
+++ b/Numerics/Numerics/Statistics/StateEstimation/ExtendedKalmanFilter.cs
@@ -1,4 +1,5 @@
-using System;
+using System;
+using CSharpNumerics.Numerics.LinearAlgebra;
using CSharpNumerics.Numerics.Objects;
namespace CSharpNumerics.Statistics.StateEstimation;
@@ -87,7 +88,8 @@ public void Update(Func h, Func jacobianH, Ma
VectorN innovation = z - h(_state);
Matrix S = H * _covariance * Ht + R;
- Matrix K = _covariance * Ht * S.Inverse();
+ // K = P·Hᵀ·S⁻¹, solved as S·Kᵀ = (P·Hᵀ)ᵀ — see KalmanFilter.Update.
+ Matrix K = S.SolveSymmetricPositiveDefinite((_covariance * Ht).Transpose()).Transpose();
_state = _state + K * innovation;
_covariance = (KalmanFilter.Identity(Dimension) - K * H) * _covariance;
diff --git a/Numerics/Numerics/Statistics/StateEstimation/KalmanFilter.cs b/Numerics/Numerics/Statistics/StateEstimation/KalmanFilter.cs
index ce660d2..cb39f5b 100644
--- a/Numerics/Numerics/Statistics/StateEstimation/KalmanFilter.cs
+++ b/Numerics/Numerics/Statistics/StateEstimation/KalmanFilter.cs
@@ -1,4 +1,5 @@
-using System;
+using System;
+using CSharpNumerics.Numerics.LinearAlgebra;
using CSharpNumerics.Numerics.Objects;
namespace CSharpNumerics.Statistics.StateEstimation;
@@ -101,7 +102,10 @@ public void Update(Matrix H, Matrix R, VectorN z)
VectorN innovation = z - H * _state; // y
Matrix S = H * _covariance * Ht + R; // innovation covariance
- Matrix K = _covariance * Ht * S.Inverse(); // Kalman gain (n×m)
+ // Kalman gain K = P·Hᵀ·S⁻¹ (n×m). Solved as S·Kᵀ = (P·Hᵀ)ᵀ rather than by forming
+ // S⁻¹: S is symmetric positive definite, so Cholesky applies, and an explicit inverse
+ // costs several times a factorize-and-solve while losing accuracy.
+ Matrix K = S.SolveSymmetricPositiveDefinite((_covariance * Ht).Transpose()).Transpose();
_state = _state + K * innovation;
diff --git a/Numerics/Numerics/Statistics/StateEstimation/KalmanSmoother.cs b/Numerics/Numerics/Statistics/StateEstimation/KalmanSmoother.cs
index 72d8403..0465ec5 100644
--- a/Numerics/Numerics/Statistics/StateEstimation/KalmanSmoother.cs
+++ b/Numerics/Numerics/Statistics/StateEstimation/KalmanSmoother.cs
@@ -1,5 +1,6 @@
-using System;
+using System;
using System.Collections.Generic;
+using CSharpNumerics.Numerics.LinearAlgebra;
using CSharpNumerics.Numerics.Objects;
namespace CSharpNumerics.Statistics.StateEstimation;
@@ -90,7 +91,8 @@ public KalmanSmootherResult Smooth(
// Update
VectorN innovation = measurements[k] - H * xPred;
Matrix S = H * pPred * Ht + R;
- Matrix K = pPred * Ht * S.Inverse();
+ // K = P_pred·Hᵀ·S⁻¹, solved as S·Kᵀ = (P_pred·Hᵀ)ᵀ.
+ Matrix K = S.SolveSymmetricPositiveDefinite((pPred * Ht).Transpose()).Transpose();
VectorN xFilt = xPred + K * innovation;
Matrix pFilt = (identity - K * H) * pPred;
@@ -112,7 +114,11 @@ public KalmanSmootherResult Smooth(
for (int k = n - 2; k >= 0; k--)
{
// Smoother gain: C = P_f[k] Fᵀ (P_pred[k+1])⁻¹
- Matrix C = filteredCovariances[k] * Ft * predictedCovariances[k + 1].Inverse();
+ // C = P_f[k]·Fᵀ·P_pred[k+1]⁻¹, solved as P_pred[k+1]·Cᵀ = (P_f[k]·Fᵀ)ᵀ.
+ // A predicted covariance is symmetric positive definite for the same reason S is.
+ Matrix C = predictedCovariances[k + 1]
+ .SolveSymmetricPositiveDefinite((filteredCovariances[k] * Ft).Transpose())
+ .Transpose();
smoothedStates[k] = filteredStates[k]
+ C * (smoothedStates[k + 1] - predictedStates[k + 1]);
diff --git a/docs/AdvancedGameEngineRoadMap.md b/docs/AdvancedGameEngineRoadMap.md
deleted file mode 100644
index 041d063..0000000
--- a/docs/AdvancedGameEngineRoadMap.md
+++ /dev/null
@@ -1,218 +0,0 @@
-# Advanced Game Engine Roadmap
-
-## Goal
-
-Extend the existing `Engines/Game` engine into a high-fidelity simulation platform targeting:
-- **Flight simulators** (6DOF aircraft, aerodynamics, atmosphere)
-- **Real-time fluid simulations** (wind, smoke, water — game-quality CFD)
-- **ML-driven game AI** (NPC behavior, adaptive difficulty, physics-informed control)
-
-The library will be distributed as Unity-compatible assets for developers building advanced simulation games.
-
-## Architecture Principle
-
-The engine **does not** replace or duplicate code in Numerics, Physics, or ML. It orchestrates them:
-
-```
-Engines.Game.Flight → Physics.FluidDynamics + Physics.Mechanics + Numerics.ODE
-Engines.Game.Fluids → Physics.FluidDynamics + Numerics.FiniteDifference
-Engines.Game.AI → ML.ReinforcementLearning + ML.NeuralNetwork
-Engines.Game (existing) → Physics.Mechanics (RigidBody, collisions, constraints)
-```
-
-Namespace: `CSharpNumerics.Engines.Game.*`
-
----
-
-## Phase 1 — Flight Dynamics Core
-
-Flight simulation foundation: 6DOF rigid-body aircraft model with aerodynamic forces.
-
-### Numerics additions (Numerics namespace)
-- [ ] Quaternion-based rotation utilities (if not already sufficient — verify QuaternionAlgebra coverage for body↔world frame transforms)
-- [ ] Frame transformation helpers (body frame ↔ world frame, wind frame ↔ body frame)
-
-### Physics additions (Physics namespace)
-- [ ] `Physics.FluidDynamics.Aerodynamics.AtmosphereModel` — ISA standard atmosphere (density, pressure, temperature vs altitude)
-- [ ] `Physics.FluidDynamics.Aerodynamics.AirfoilModel` — lift/drag coefficient tables (AoA-based lookup with interpolation, flat plate, NACA symmetric)
-- [ ] `Physics.FluidDynamics.Aerodynamics.ControlSurface` — deflection → ΔCl/ΔCd model (elevator, aileron, rudder)
-- [ ] `Physics.FluidDynamics.Aerodynamics.PropulsionModel` — thrust as function of throttle, altitude, speed (jet + propeller)
-
-### Engine additions (Engines.Game namespace)
-- [ ] `Engines.Game.Flight.AircraftState` — 12-state vector (position, velocity, quaternion attitude, angular rates)
-- [ ] `Engines.Game.Flight.FlightDynamicsEngine` — implements `ISimulationEngine`; integrates 6DOF equations using RK4Stepper from Numerics
-- [ ] `Engines.Game.Flight.ControlInput` — throttle, pitch, roll, yaw, flaps, gear
-- [ ] `Engines.Game.Flight.AircraftConfig` — wing area, span, mass, CG position, engine count, control surface geometry
-
-### Tests
-- [ ] Straight-and-level flight (trim condition)
-- [ ] Stall behavior (AoA > critical → Cl drops)
-- [ ] Climb/descent rates match expected performance
-- [ ] Roll/pitch response to control inputs
-
----
-
-## Phase 2 — Real-Time Fluid Simulation for Games
-
-Adapt existing Navier-Stokes solvers for game-quality real-time performance. Focus on visual fidelity at interactive framerates rather than engineering accuracy.
-
-### Physics additions (Physics namespace)
-- [ ] `Physics.FluidDynamics.Turbulence.TurbulenceModel` — simple k-epsilon or Smagorinsky SGS for visual turbulence
-- [ ] `Physics.FluidDynamics.Buoyancy.BuoyancyForce` — hot gas rises, cold sinks (for smoke/fire)
-
-### Engine additions (Engines.Game namespace)
-- [ ] `Engines.Game.Fluids.GameFluidSolver2D` — stripped-down 2D N-S optimized for speed (fewer Poisson iterations, coarser grid, LOD)
-- [ ] `Engines.Game.Fluids.GameFluidSolver3D` — 3D variant with configurable resolution (32³–128³)
-- [ ] `Engines.Game.Fluids.FluidConfig` — grid size, viscosity, timestep, boundary mode, quality preset (Low/Medium/High)
-- [ ] `Engines.Game.Fluids.FluidBodyCoupling` — two-way: rigid bodies experience drag/lift from fluid; bodies displace fluid
-- [ ] `Engines.Game.Fluids.VorticityConfinement` — add back vorticity lost to numerical diffusion (sharper smoke curls)
-- [ ] `Engines.Game.Fluids.FluidEmitter` — inject velocity/density at a point or region (exhaust, wind sources, explosions)
-- [ ] `Engines.Game.Fluids.FluidObstacle` — static/dynamic geometry that blocks flow (buildings, terrain, aircraft)
-
-### Integration with Flight
-- [ ] Aircraft samples wind field at its position → adds to aerodynamic velocity
-- [ ] Engine exhaust feeds back into fluid as emitter
-- [ ] Wake turbulence behind aircraft (simplified vortex model)
-
-### Tests
-- [ ] Smoke rising around obstacle (visual validation: vortex shedding)
-- [ ] Fluid-body coupling: sphere falls slower in fluid vs vacuum
-- [ ] Emitter produces expanding plume
-- [ ] Performance: 64³ grid runs at ≥30 steps/sec on single thread
-
----
-
-## Phase 3 — ML-Driven Game AI
-
-Wire existing RL infrastructure into the game engine for intelligent NPC behavior and adaptive systems.
-
-### ML additions (ML namespace)
-- [ ] `ML.ReinforcementLearning.Environments.FlightEnv` — observation: aircraft state (12D), action: control surfaces (4D continuous), reward: waypoint tracking + fuel efficiency
-- [ ] `ML.ReinforcementLearning.Environments.DogfightEnv` — multi-agent pursuit-evasion with flight dynamics
-- [ ] `ML.ReinforcementLearning.Environments.FluidNavigationEnv` — agent navigates through wind/current field
-
-### Engine additions (Engines.Game namespace)
-- [ ] `Engines.Game.AI.GameAIAgent` — wraps a trained RL policy; takes game state observations, returns actions per tick
-- [ ] `Engines.Game.AI.AITrainer` — offline training loop: runs FlightDynamicsEngine headless, trains PPO/DDPG policy
-- [ ] `Engines.Game.AI.BehaviorTree` — lightweight behavior tree executor for hybrid AI (ML decisions + scripted fallbacks)
-- [ ] `Engines.Game.AI.FormationController` — multi-agent coordination (wingmen following leader using RL + formation rules)
-- [ ] `Engines.Game.AI.AdaptiveDifficulty` — monitors player performance, adjusts AI aggressiveness via policy parameters
-
-### Tests
-- [ ] Trained agent can maintain level flight (reward converges)
-- [ ] Dogfight agent pursues target with reasonable intercept geometry
-- [ ] Formation controller maintains spacing under wind perturbation
-- [ ] Adaptive difficulty reduces AI skill when player struggles
-
----
-
-## Phase 4 — Advanced Physics Integration
-
-Deeper physics features needed for high-fidelity simulation games.
-
-### Physics additions (Physics namespace)
-- [ ] `Physics.Mechanics.SoftBody.DeformableMesh` — mass-spring network on triangulated mesh (uses CoupledOscillators)
-- [ ] `Physics.Mechanics.SoftBody.ClothSimulation` — constrained particle system with self-collision
-- [ ] `Physics.FluidDynamics.SPH.SPHSolver` — Smoothed Particle Hydrodynamics for water splashes, liquid in containers
-- [ ] `Physics.FluidDynamics.FreeSurface.VOFTracker` — Volume-of-Fluid for water surface tracking
-
-### Engine additions (Engines.Game namespace)
-- [ ] `Engines.Game.Particles.ParticleSystem` — emission, lifetime, forces (gravity, drag, wind from fluid field), collision with world
-- [ ] `Engines.Game.Particles.ParticleEmitter` — rate, cone angle, initial velocity, randomness
-- [ ] `Engines.Game.Terrain.TerrainCollider` — heightmap-based collision for ground interaction
-- [ ] `Engines.Game.Terrain.WindOverTerrain` — terrain deflects wind field (simple boundary layer model)
-- [ ] Continuous collision detection (CCD) — swept sphere for fast projectiles
-- [ ] Spatial partitioning upgrade — BVH or octree for large worlds
-
-### Tests
-- [ ] Cloth draped over sphere (visual: no penetration, natural drape)
-- [ ] SPH water splashes when object dropped
-- [ ] Particle system affected by fluid wind field
-- [ ] Terrain collision for aircraft landing gear
-
----
-
-## Phase 5 — Unity Integration Layer
-
-Package the engine for consumption as a Unity asset. This layer is a thin adapter — all logic stays in the core library.
-
-### Engine additions (Engines.Game namespace)
-- [ ] `Engines.Game.Unity.UnityAdapter` — converts between Unity Vector3/Quaternion and CSharpNumerics VectorN/Matrix
-- [ ] `Engines.Game.Unity.PhysicsSync` — synchronizes PhysicsWorld state with Unity Transform components (position, rotation)
-- [ ] `Engines.Game.Unity.FluidRenderer` — provides density/velocity textures from GameFluidSolver3D for Unity VFX Graph or shader consumption
-- [ ] `Engines.Game.Unity.FlightController` — MonoBehaviour-style interface: exposes ControlInput, reads AircraftState, drives Unity Transform
-- [ ] `Engines.Game.Unity.AIBridge` — feeds Unity game state to GameAIAgent, applies returned actions
-
-### Documentation & Examples
-- [ ] Unity sample project: flight simulator scene with HUD
-- [ ] Unity sample: real-time smoke simulation with fluid-body interaction
-- [ ] Unity sample: ML-trained dogfight AI opponents
-- [ ] API reference documentation for all public types
-- [ ] Update Engines/README.md with Game Engine section
-
-### Tests
-- [ ] Adapter round-trips preserve values (Vector3 → VectorN → Vector3)
-- [ ] PhysicsSync maintains frame coherence over 10,000 steps
-- [ ] FluidRenderer produces valid texture data at 60fps
-
----
-
-## Phase 6 — Performance & Polish
-
-Optimize for production game use.
-
-- [ ] SIMD acceleration for vector math hot paths (Vector128/256 where available)
-- [ ] Span / stackalloc for zero-allocation fluid solver inner loops
-- [ ] Multi-threaded broad phase collision (partition spatial grid across threads)
-- [ ] Fluid solver thread: run on background thread, game samples latest snapshot
-- [ ] Memory pooling for particles, rigid bodies, constraint arrays
-- [ ] Profiling benchmarks: measure ms/frame for key scenarios
-- [ ] LOD system for fluid: coarse grid far from camera, fine near
-- [ ] Deterministic replay: serialize full state for network sync / replays
-
----
-
-## Summary — What Goes Where
-
-| Feature | Namespace | Rationale |
-|---------|-----------|-----------|
-| Atmosphere model | `Physics.FluidDynamics.Aerodynamics` | Physical model, no game logic |
-| Airfoil lift/drag | `Physics.FluidDynamics.Aerodynamics` | Physical model |
-| Turbulence model | `Physics.FluidDynamics.Turbulence` | Physical model |
-| SPH solver | `Physics.FluidDynamics.SPH` | Physical algorithm |
-| Soft body / cloth | `Physics.Mechanics.SoftBody` | Physical model |
-| Quaternion math | `Numerics.Objects` | Pure math |
-| Flight dynamics engine | `Engines.Game.Flight` | Orchestration |
-| Game fluid solver | `Engines.Game.Fluids` | Game-optimized orchestration |
-| Game AI agent | `Engines.Game.AI` | Game-specific ML application |
-| Particle system | `Engines.Game.Particles` | Game feature |
-| Terrain interaction | `Engines.Game.Terrain` | Game feature |
-| Unity adapter | `Engines.Game.Unity` | Platform integration |
-| RL flight environment | `ML.ReinforcementLearning.Environments` | Training environment (ML domain) |
-
----
-
-## Decision: New Engine or Extend Existing?
-
-**Extend the existing `Engines/Game` engine.** Rationale:
-
-1. The current PhysicsWorld + constraint solver + collision pipeline is solid and tested
-2. Flight, fluids, and AI are **additive sub-systems** — they don't replace rigid-body physics
-3. The `ISimulationEngine` + `SimulationClock` + `EventBus` infrastructure already supports multiple coordinated engines
-4. Creating a separate engine would duplicate the rigid-body foundation and violate DRY
-5. Unity integration is just an adapter layer — it doesn't require engine-level changes
-
-The extended structure:
-```
-Engines/Game/
-├── PhysicsWorld.cs (existing — rigid body orchestrator)
-├── Flight/ (NEW — 6DOF flight dynamics)
-├── Fluids/ (NEW — real-time game fluid sim)
-├── AI/ (NEW — ML-driven game intelligence)
-├── Particles/ (NEW — particle effects)
-├── Terrain/ (NEW — terrain interaction)
-├── Unity/ (NEW — Unity asset bridge)
-├── Constraints/ (existing)
-├── BroadPhase/ (existing)
-└── Objects/ (existing)
-```
diff --git a/docs/ExoplanetEngineRoadMap.md b/docs/ExoplanetEngineRoadMap.md
deleted file mode 100644
index 94eef04..0000000
--- a/docs/ExoplanetEngineRoadMap.md
+++ /dev/null
@@ -1,244 +0,0 @@
-# Exoplanet Engine Roadmap
-
-## Mål
-
-Bygga `CSharpNumerics.Engines.Exoplanet` — en analysmotor för exoplanet-transitdetektion som kombinerar befintliga moduler:
-
-- **`CSharpNumerics.Physics.Optics`** — limbmörkning, transitgeometri
-- **`CSharpNumerics.Statistics.Fitting`** — transitmodell-anpassning (Mandel-Agol, trapetsmodell)
-- **`CSharpNumerics.Statistics.TimeSeriesAnalysis`** — BLS, Lomb-Scargle, avtrending, fasveckning
-- **`CSharpNumerics.ML.Sequence`** — CNN1D/LSTM/BiLSTM-klassificering av ljuskurvefönster
-
-**Slutmål:** Träna på märkt (labeled) transitdata → exportera tränad modell → prediktioner via webbtjänst.
-
-API-integrationer mot NASA Exoplanet Archive, MAST/TESS, etc. sker **utanför** CSharpNumerics, men interna datamodeller ska matcha standardformaten (Kepler/TESS FITS-kolumner, KOI-tabeller).
-
----
-
-## Nulägesanalys — Befintliga resurser
-
-| Komponent | Status | Namespace / Klass |
-|-----------|--------|-------------------|
-| BLS-periodanalys | ✅ | `Statistics.TimeSeriesAnalysis.BoxFittingLeastSquares` |
-| Lomb-Scargle-periodogram | ✅ | `Statistics.TimeSeriesAnalysis.LombScarglePeriodogram` |
-| Fasveckning & binning | ✅ | `Statistics.TimeSeriesAnalysis.PhaseFolding` |
-| Avtrending (polynom, SG, median) | ✅ | `Statistics.TimeSeriesAnalysis.TimeSeriesDetrending` |
-| Toppdetektering | ✅ | `Statistics.TimeSeriesAnalysis.PeakFitting` |
-| Icke-linjär minsta-kvadrat (LM) | ✅ | `Statistics.Fitting.NonlinearLeastSquaresFitter` |
-| Robust anpassning (IRLS) | ✅ | `Statistics.Fitting.RobustFitter` |
-| Bootstrap-konfidensintervall | ✅ | `Statistics.Fitting.ParameterEstimation` |
-| Modellval (AIC/BIC) | ✅ | `Statistics.Fitting.GoodnessOfFit` |
-| CNN1DClassifier | ✅ | `ML.Sequence.Models.Classification.CNN1DClassifier` |
-| LSTMClassifier | ✅ | `ML.Sequence.Models.Classification.LSTMClassifier` |
-| BiLSTMClassifier | ✅ | `ML.Sequence.Models.Classification.BiLSTMClassifier` |
-| SupervisedExperiment (grid search, CV) | ✅ | `ML.Experiment.SupervisedExperiment` |
-| SequenceDataHelper.CreateWindows | ✅ | `ML.Sequence.SequenceDataHelper` |
-| TimeSeries-datastruktur | ✅ | `Statistics.Data.TimeSeries` |
-| Ray-tracing, optiska ytor, sensor | ✅ | `Physics.Optics.*` |
-| Astronomi (avstånd-konverteringar) | ⚠️ Minimal | `Physics.Astro.AstronomyExtensions` |
-| Engine-infrastruktur (ISimulationEngine, EventBus, Clock) | ✅ | `Engines.Common.*` |
-| Exoplanet-specifik transitfysik | ❌ | — |
-| Ljuskurve-datamodell (TESS/Kepler-kompatibel) | ❌ | — |
-| Feature-extraktion för transiter | ❌ | — |
-| Inference-pipeline (träna → exportera → prediktera) | ❌ | — |
-
----
-
-## Externt dataformat — Referens
-
-Integrationer sker utanför CSharpNumerics, men interna typer ska mappas **1:1** mot standardformat. Nedan sammanfattas de viktigaste kolumnerna.
-
-### Kepler/TESS ljuskurva (tidsserieformat)
-
-| Kolumn | Typ | Enhet | Beskrivning |
-|--------|-----|-------|-------------|
-| `TIME` | `double` | BJD − 2457000 (TESS) / BJD − 2454833 (Kepler) | Barycentrisk Julian Date |
-| `SAP_FLUX` | `double` | e⁻/s | Simple Aperture Photometry |
-| `PDCSAP_FLUX` | `double` | e⁻/s | Pre-search Data Conditioned flux (systematikrensad) |
-| `SAP_FLUX_ERR` | `double` | e⁻/s | Osäkerhet SAP |
-| `PDCSAP_FLUX_ERR` | `double` | e⁻/s | Osäkerhet PDCSAP |
-| `QUALITY` | `int` | flaggmask | Datakvalitetsflaggor (0 = ok) |
-| `CADENCENO` | `int` | — | Sekvensnummer |
-
-### KOI / TOI transit-parametrar
-
-| Parameter | Kolumn (Kepler) | Kolumn (TESS TOI) | Enhet |
-|-----------|-----------------|---------------------|-------|
-| Orbital period | `koi_period` | `pl_orbper` | dagar |
-| Transit epoch | `koi_time0bk` | `pl_tranmid` | BJD |
-| Transit depth | `koi_depth` | `pl_trandep` | ppm |
-| Transit duration | `koi_duration` | `pl_trandur` | timmar |
-| Planet/star radius-ratio | `koi_ror` | `pl_ratror` | — |
-| Impact parameter | `koi_impact` | `pl_imppar` | — |
-| Disposition | `koi_disposition` | `tfopwg_disp` | CANDIDATE / CONFIRMED / FALSE POSITIVE |
-| Disposition score | `koi_score` | — | 0–1 |
-| Stellar Teff | `koi_steff` | `st_teff` | K |
-| Stellar radius | `koi_srad` | `st_rad` | R☉ |
-| Stellar logg | `koi_slogg` | `st_logg` | cgs |
-
----
-
-## Phase 1 — Datamodell & ljuskurvetyper
-
-Definiera de centrala datatyperna i `Engines/Exoplanet/` som matchar externt format. Inga externa beroenden — bara POCOs och konverteringslogik mot befintligt `TimeSeries`.
-
-- [x] Skapa `Engines/Exoplanet/` mappstruktur med `Data/`, `Pipeline/`, `Features/`, `Interfaces/`, `Enums/`
-- [x] `LightCurve` — klass med `double[] Time`, `double[] Flux`, `double[] FluxError`, `int[] QualityFlags`, `LightCurveMetadata Metadata`; factory-metod `FromTimeSeries(TimeSeries ts, string timeCol, string fluxCol, string errCol)`
-- [x] `LightCurveMetadata` — `string TargetId`, `string Mission` (Kepler/TESS/K2), `double TimeOffset` (BJD-offset), `CadenceType Cadence` (Short/Long)
-- [x] `TransitParameters` — record/klass: `double Period`, `double Epoch`, `double Depth`, `double Duration`, `double RadiusRatio`, `double ImpactParameter`, `double IngressDuration`
-- [x] `TransitCandidate` — `TransitParameters Parameters`, `double Score`, `TransitDisposition Disposition`, `LightCurve PhaseFoldedCurve`, `TransitFeatureSet Features`
-- [x] `StellarProperties` — `double EffectiveTemp`, `double Radius`, `double Mass`, `double SurfaceGravity`, `double Metallicity`, `SpectralType Type`
-- [x] `TransitDisposition` (enum) — `Unknown`, `Candidate`, `Confirmed`, `FalsePositive`, `AstrophysicalFalsePositive`
-- [x] `CadenceType` (enum) — `Short`, `Long`, `Fast` (TESS 20s)
-- [x] `LightCurveSanitizer` — statisk klass: `RemoveBadQuality(LightCurve, qualityMask)`, `RemoveOutliers(LightCurve, sigmaThreshold)`, `FillGaps(LightCurve, maxGapSize)`, `NormalizeFlux(LightCurve)` → returnerar ny `LightCurve`
-- [x] Enhetstester för samtliga datatyper och `LightCurveSanitizer`
-
----
-
-## Phase 2 — Transitfysik (Physics.Astro)
-
-Utöka `CSharpNumerics.Physics.Astro` med fysikaliska modeller för transit-geometri och limbmörkning. Dessa är domänoberoende fysik och tillhör Physics-sektionen.
-
-- [x] `TransitGeometry` — statisk klass i `Physics/Astro/`: `ImpactParameter(a, inclination, Rstar)`, `TransitProbability(a, Rstar, Rplanet)`, `TransitDuration(period, a, Rstar, Rplanet, inclination)`, `IngressDuration(...)`, `ContactTimes(...)` (T1–T4)
-- [x] `LimbDarkening` — statisk klass i `Physics/Astro/`: `Linear(mu, u1)`, `Quadratic(mu, u1, u2)`, `NonlinearFourParam(mu, c1..c4)`, `IntensityProfile(model, u_params, muArray)` — returnerar intensitetsprofil
-- [x] `TransitModel` — klass i `Physics/Astro/`: implementerar Mandel & Agol (2002) analytisk transitljuskurva. `Evaluate(double[] times, TransitParameters p, LimbDarkeningModel ld, StellarProperties star)` → `double[] modelFlux`. Stödjer cirkulär bana.
-- [x] `LimbDarkeningModel` (enum) i `Physics/Astro/` — `Uniform`, `Linear`, `Quadratic`, `NonlinearFourParam`
-- [x] `KeplerOrbit` — statisk klass i `Physics/Astro/`: `TrueAnomaly(meanAnomaly, eccentricity)` (Kepler-ekvationen, Newton-iteration), `SemiMajorAxis(period, stellarMass)` (Keplers tredje lag), `OrbitalVelocity(a, period)`
-- [x] Enhetstester: verifiera transitmodell mot kända Kepler-transitkurvor (analytiska testfall), limbmörkningsprofiler, Kepler-ekvationen
-
----
-
-## Phase 3 — Transitdetektionspipeline
-
-Kombinera befintliga `TimeSeriesAnalysis`-moduler till en sammanhållen detektionspipeline i `Engines/Exoplanet/Pipeline/`.
-
-- [x] `LightCurvePreprocessor` — orkestrerar: `Detrend(LightCurve, method)` → `RemoveOutliers()` → `Normalize()`. Delegerar till `TimeSeriesDetrending` och `LightCurveSanitizer`
-- [x] `PeriodSearcher` — wrapper kring BLS och Lomb-Scargle: `Search(LightCurve, minPeriod, maxPeriod, method)` → `PeriodSearchResult` (bästa period, FAP, spektrum). Stödjer `PeriodSearchMethod`-enum (`BLS`, `LombScargle`, `Both`)
-- [x] `TransitFitter` — anpassar `TransitModel` (Phase 2) till fasveckt ljuskurva via `NonlinearLeastSquaresFitter`. `Fit(LightCurve, trialPeriod, trialEpoch)` → `TransitFitResult` med `TransitParameters`, osäkerheter, χ²/BIC
-- [x] `TransitValidator` — utvärderar kandidater: kontrollerar SNR (transit-djup/brus), udda/jämn transit-djup-test, sekundär eklipsanalys, centrumfotoanalys. `Validate(TransitCandidate, LightCurve)` → `ValidationResult` med `IsValid`, `Warnings[]`, `Score`
-- [x] `TransitDetectionPipeline` — main orchestrator: `Detect(LightCurve, config)` → `TransitCandidate[]`. Steg: preprocessor → period search → phase fold → fit → validate → upprepa (multi-planet med iterativ subtraktion)
-- [x] `TransitDetectionConfig` — konfigurationsobjekt: `MinPeriodDays`, `MaxPeriodDays`, `MinTransitDepthPpm`, `SnrThreshold`, `MaxPlanets`, `DetrendingMethod`, `PeriodSearchMethod`
-- [x] Enhetstester: syntetisk ljuskurva med 1–3 injicerade transiter → pipeline ska återfinna samtliga med korrekt period (±1%)
-
----
-
-## Phase 4 — Feature-extraktion
-
-Extrahera transitspecifika features för ML-klassificering. Placeras i `Engines/Exoplanet/Features/`.
-
-- [x] `TransitFeatureExtractor` — extraherar features från en `TransitCandidate` + rå `LightCurve`: transit-djup, duration, ingress/egress-ratio, periodik-styrka (BLS SNR), udda/jämn djupskillnad, sekundäreklips-djup, ljuskurve-scatter, fasvecknings-χ², limbmörkningskoefficienter (anpassade), transitform-metric (V-shape vs U-shape)
-- [x] `TransitFeatureSet` — resultatobjekt: `Dictionary Features`, med standardnamn: `Depth`, `Duration`, `Period`, `SnrBls`, `OddEvenRatio`, `VShapeMetric`, `SecondaryDepth`, `ScatterInTransit`, `ScatterOutTransit`, `LimbDarkeningU1`, `LimbDarkeningU2`, `IngressEgressRatio`
-- [x] `WindowedFeatureExtractor` — skapar `(Matrix X, VectorN y)` för ML-träning: kombinerar raw fasveckt ljuskurva (tidsfönster) med extraherade features som extra kolumner. Kompatibel med `SequenceDataHelper.CreateWindows`
-- [x] Enhetstester: verifiera extraherade features mot kända transiter
-
----
-
-## Phase 5 — ML-träning & inference
-
-Bygg tränings- och inference-workflow i `Engines/Exoplanet/Pipeline/`. Målet: träna en modell på labeled data → serialisera → ladda för prediktion.
-
-- [x] `TransitClassifierTrainer` — façade som sammanfogar feature-extraktion + `SupervisedExperiment`. `Train(LightCurve[] curves, TransitDisposition[] labels, TrainerConfig config)` → `TrainedTransitModel`. Grid-search över CNN1D, LSTM, BiLSTM med konfigurerbar hyperparameterrymd
-- [x] `TrainerConfig` — konfiguration: `int WindowSize`, `int Stride`, `ModelType[] CandidateModels`, `CrossValidatorConfig Cv`, `int Epochs`, `double ValidationSplit`
-- [x] `TrainedTransitModel` — wrapper: innehåller den tränade `IClassificationModel`, `TransitDetectionConfig`, feature-scaler-parametrar, träningsmetrik (accuracy, precision, recall, F1, confusion matrix). Stödjer `Predict(LightCurve)` → `TransitPrediction[]`
-- [x] `TransitPrediction` — `TransitCandidate Candidate`, `double Probability`, `TransitDisposition PredictedDisposition`
-- [x] `ModelSerializer` — serialisering/deserialisering av `TrainedTransitModel` till/från byte-array eller fil. Använder `FieldSerializer` från `Engines.Common` som bas. Format: JSON-metadata + binär viktdata
-- [x] `TransitInferencePipeline` — stateless inference: `LoadModel(byte[] serialized)` → `Predict(LightCurve)` → `TransitPrediction[]`. Designad för att anropas från webbtjänst. Trådsäker (inga mutable state)
-- [x] Enhetstester: full roundtrip — träna på syntetisk data → serialisera → deserialisera → prediktera → verifiera accuracy, precision, recall
-
----
-
-## Phase 6 — Statistiska tillägg (Statistics)
-
-Nya statistiska metoder som behövs för transitanalys men är generellt användbara. Placeras i respektive Statistics-namespace.
-
-- [x] `SigmaClipping` i `Statistics/Robust/` (`CSharpNumerics.Statistics.Robust`) — iterativ sigma-clipping: `Clip(double[] data, double sigmaLow, double sigmaHigh, int maxIter)` → `ClipResult` (mask, mean, std). Används av `LightCurveSanitizer.RemoveOutliers`
-- [x] `FalseAlarmProbability` i `Statistics/TimeSeriesAnalysis/` — FAP-beräkning för periodogram-toppar: analytisk (Baluev 2008) och bootstrap-baserad. `AnalyticalFAP(peakPower, numFrequencies, numDataPoints)`, `BootstrapFAP(times, values, peakPower, nBootstrap)`
-- [x] `SlidingWindowStatistics` i `Statistics/Robust/` (`CSharpNumerics.Statistics.Robust`) — löpande medelvärde, standardavvikelse, MAD (median absolute deviation) med konfigurbar fönsterstorlek. Används för scatter-beräkning i transit-features
-- [x] Enhetstester för samtliga statistiska tillägg
-
----
-
-## Phase 7 — Engine-integration & demo
-
-Integrera allt till `ExoplanetEngine : ISimulationEngine` och bygg en end-to-end demo.
-
-- [x] `ExoplanetEngine` — implementerar `ISimulationEngine`. `Init()` laddar konfiguration och eventuellt tränad modell. `Step(dt)` processar nästa ljuskurva-segment. Publicerar `TransitDetectedEvent` via `EventBus` vid fynd
-- [x] `TransitDetectedEvent` — eventtyp: `TransitCandidate Candidate`, `double Timestamp`
-- [x] `ExoplanetEngineConfig` — konfiguration: `TransitDetectionConfig Detection`, `TrainerConfig Training`, `string ModelPath` (för förtränad modell)
-- [x] End-to-end integrationstest: syntetisera TESS-liknande ljuskurvor (realistisk kadens, brus, systematics) → kör full pipeline (avtrend → BLS → fit → ML-klassificering) → verifiera att korrekt antal transiterande planeter detekteras
-- [x] Demo-test som visar hela workflow: data-inläsning → träning → inference → resultatutskrift. Dokumenterar det format som en webbtjänst behöver skicka/ta emot
-- [x] Uppdatera `Engines/README.md` med Exoplanet Engine-dokumentation: arkitekturdiagram, API-referens, kodexempel, dataformatbeskrivning
-
----
-
-## Sammanfattning — Berörda namespaces
-
-| Namespace | Typ av ändring |
-|-----------|---------------|
-| `Engines.Exoplanet` | **Ny** — motor, pipeline, features, datamodell |
-| `Engines.Exoplanet.Data` | **Ny** — `LightCurve`, `TransitCandidate`, `StellarProperties`, etc. |
-| `Engines.Exoplanet.Pipeline` | **Ny** — detektion, träning, inference |
-| `Engines.Exoplanet.Features` | **Ny** — feature-extraktion |
-| `Physics.Astro` | **Utökad** — `TransitGeometry`, `LimbDarkening`, `TransitModel`, `KeplerOrbit`, `LimbDarkeningModel` (enum) |
-| `Statistics.Robust` | **Ny** — `SigmaClipping` |
-| `Statistics` | **Ny** — `SlidingWindowStatistics` |
-| `Statistics.TimeSeriesAnalysis` | **Utökad** — `FalseAlarmProbability` |
-| `ML.Sequence` | Befintlig — används som är |
-| `Statistics.Fitting` | Befintlig — används som är |
-| `Engines.Common` | Befintlig — `ISimulationEngine`, `EventBus`, `FieldSerializer` |
-
----
-
-## Dataflöde — Arkitekturskiss
-
-```
- ┌─────────────────────────────────────────────┐
- │ EXTERN WEBBTJÄNST │
- │ (hämtar data från NASA/MAST/TESS API) │
- └──────────────────┬──────────────────────────┘
- │ LightCurve (JSON/CSV)
- ▼
- ┌─────────────────────────────────────────────┐
- │ Engines.Exoplanet │
- │ │
- │ ┌─────────────┐ ┌────────────────────┐ │
- │ │ LightCurve │──▶│ LightCurvePreproc │ │
- │ │ Sanitizer │ │ (Detrend/Normalize)│ │
- │ └─────────────┘ └────────┬───────────┘ │
- │ │ │
- │ ┌────────────▼────────────┐ │
- │ │ PeriodSearcher │ │
- │ │ (BLS / Lomb-Scargle) │ │
- │ └────────────┬────────────┘ │
- │ │ │
- │ ┌────────────▼────────────┐ │
- │ │ TransitFitter │ │
- │ │ (Mandel-Agol via LM) │ │
- │ └────────────┬────────────┘ │
- │ │ │
- │ ┌────────────────┐ │ │
- │ │FeatureExtractor│◀───┘ │
- │ └───────┬────────┘ │
- │ │ │
- │ ┌────────────▼──────────────────────┐ │
- │ │ ML Inference (CNN1D/LSTM/BiLSTM) │ │
- │ │ via TrainedTransitModel │ │
- │ └────────────┬──────────────────────┘ │
- │ │ │
- │ ┌────────────▼────────────┐ │
- │ │ TransitValidator │ │
- │ │ (SNR, odd/even, etc.) │ │
- │ └────────────┬────────────┘ │
- │ │ │
- │ ▼ │
- │ TransitPrediction[] │
- │ (→ JSON till webbtjänst) │
- └─────────────────────────────────────────────┘
-
-Beroenden:
- Engines.Exoplanet → Physics.Astro (TransitModel, LimbDarkening, KeplerOrbit)
- Engines.Exoplanet → Statistics.Fitting (NonlinearLeastSquaresFitter)
- Engines.Exoplanet → Statistics.TimeSeriesAnalysis (BLS, LombScargle, Detrending)
- Engines.Exoplanet → ML.Sequence (CNN1DClassifier, LSTMClassifier, BiLSTMClassifier)
- Engines.Exoplanet → Engines.Common (ISimulationEngine, EventBus, FieldSerializer)
-```
diff --git a/docs/Multiphysics-Roadmap.md b/docs/Multiphysics-Roadmap.md
deleted file mode 100644
index 9fb7cd5..0000000
--- a/docs/Multiphysics-Roadmap.md
+++ /dev/null
@@ -1,140 +0,0 @@
-# Multiphysics Engine Roadmap
-
-## Overview
-
-A new `Engines/Multiphysics/` simulation engine enabling four student-friendly PDE simulations on simple geometries, designed for visualization in Unity or web platforms. Built on existing finite difference infrastructure with new engineering materials, Unity/web export, and optional ML/Monte Carlo integration following the GIS engine pattern.
-
-**Target complexity:** ~20–25 files, ~3.5–5K LOC (comparable to Game/Quantum engines).
-
-### LBM Assessment
-
-Lattice Boltzmann Methods were evaluated and **rejected**. FD (Grid2D, GridOperators, 4 time steppers) already covers all four simulation types. LBM would require ~1–2K LOC of new infrastructure (D2Q9 distributions, BGK collision, streaming) duplicating existing FD capability. LBM strengths (parallelism, complex boundaries) are unnecessary for simple student geometries.
-
-### Simulation Types
-
-| Simulation | PDE | Method | Dimension |
-|------------|-----|--------|-----------|
-| Heat plate | ∂T/∂t = α∇²T | FD (Grid2D + Laplacian2D + ITimeStepper) | 2D |
-| Pipe flow | Navier-Stokes (simplified) | FD (1D transient Hagen-Poiseuille) | 1D (2D stretch) |
-| Electric field | ∇²φ = −ρ/ε | FD (Poisson solver + Gradient2D) | 2D |
-| Beam stress | EIu⁗ = q | 1D FEM (Euler-Bernoulli) | 1D |
-
----
-
-## Phase 1 — Foundations & Shared Infrastructure
-
-- [x] Add iterative Poisson solver (Gauss-Seidel) to `GridOperators` for ∇²φ = f on Grid2D with Dirichlet BC — `GridOperators.SolvePoisson2D()`
-- [x] Add `SolidExtensions` to `Physics/` — Hooke's law, second moment of area, Euler-Bernoulli beam equation, analytical deflections
-- [x] Add minimal 1D FEM primitives in `Numerics/Numerics/Numerics/FiniteElement/`
- - [x] `IElement1D` interface — shape functions, local stiffness, local load
- - [x] `BarElement` (2-node, linear) — axial stress/strain
- - [x] `BeamElement` (2-node, Hermite cubic) — Euler-Bernoulli bending
- - [x] `Assembler1D` — local→global stiffness matrix assembly with Gaussian elimination solver
- - [x] `Mesh1D` — 1D mesh from interval subdivision (nodes + elements)
-- [x] Add `EngineeringMaterial` to `Physics/Materials/Engineering/`
- - [x] `EngineeringMaterial` immutable struct: ThermalConductivity, SpecificHeat, Density, DynamicViscosity, ElectricPermittivity, YoungsModulus, PoissonsRatio
- - [x] `EngineeringLibrary` with common materials: steel, aluminum, copper, water, air, concrete, glass
- - [ ] Extend `Materials` factory with `Materials.Engineering("Steel")` — add `EngineeringMaterial?` property to `MaterialDescriptor`
-
-## Phase 2 — Core Engine & Individual Solvers
-
-- [x] Create engine scaffold in `Engines/Multiphysics/` with **fluent API** (*pattern from GIS `RiskScenario`*)
- - [x] `SimulationType` static entry point — `SimulationType.Create(MultiphysicsType)`
- - [x] `SimulationBuilder` — fluent builder: `.WithMaterial()` → `.WithGeometry()` → `.WithBoundary()` → `.WithInitialCondition()` → `.AddSource()` → `.Solve()` / `.Run()`
- - [x] `MultiphysicsType` enum: HeatPlate, PipeFlow, ElectricField, BeamStress
- - [x] `BeamSupport` enum: Cantilever, SimplySupported, FixedFixed
- - [x] `SimulationResult` — unified result type with 2D fields, 1D arrays, beam-specific outputs, E-field vectors
- - [x] `IMultiphysicsSolver` — internal solver interface
-- [x] `HeatPlateSolver` — 2D heat equation ∂T/∂t = α∇²T + source (Grid2D + Laplacian2D + forward Euler)
- - [x] Material-driven α = k/(ρ·cp) from `EngineeringMaterial`
- - [x] BCs: fixed temperature (Dirichlet edges)
- - [x] Initial conditions: uniform or function-based
-- [x] `PipeFlowSolver` — 1D transient Hagen-Poiseuille via FD (cylindrical Laplacian)
- - [x] Material viscosity ν from `EngineeringMaterial.KinematicViscosity`
- - [x] Symmetry BC at r=0, no-slip at r=R
- - [ ] Optional stretch: 2D lid-driven cavity (projection method + Poisson solver)
-- [x] `ElectricFieldSolver` — 2D Poisson/Laplace ∇²φ = −ρ/ε
- - [x] Poisson solver from Phase 1, E = −∇φ via `GridOperators.Gradient2D`
- - [x] Conductor as BC (φ = V on edge cells), charge distributions as source
- - [x] Material permittivity from `EngineeringMaterial`
-- [x] `BeamStressSolver` — 1D Euler-Bernoulli via analytical `SolidExtensions` formulas
- - [x] Support: cantilever, simply supported, fixed-fixed
- - [x] Loads: point loads, distributed loads
- - [x] Output: deflection curve, bending moment, shear force, stress distribution
- - [x] EI from `EngineeringMaterial.YoungsModulus` + cross-section (rectangular, circular, or custom I)
-
-## Phase 3 — Export & Visualization Support
-
-- [x] `MultiphysicsBinaryExporter` — Unity-compatible binary format (MPHY magic + header + float layers)
- - [x] 2D timeline export: per-step float[Nx*Ny] layers
- - [x] 1D beam export: positions + 4 layers (deflection, moment, shear, stress)
- - [x] Reader: `ReadHeader()` and `Read()` for round-trip loading
-- [x] `MultiphysicsJsonExporter` — JSON format for web visualization
- - [x] `SimulationResult` direct export (all 4 simulation types)
- - [x] `SimulationTimeline` export with per-step snapshots
- - [x] `FieldSnapshot` single-frame export
- - [x] `BeamSnapshot` export with all curves
- - [x] `ExportMetadata` support (simulation, unit, density)
- - [x] `Save()` to file for all types
-- [x] Snapshot system
- - [x] `FieldSnapshot` — time-stamped 2D scalar field with flat indexing, Min/Max, ToArray round-trip
- - [x] `BeamSnapshot` — immutable 1D arrays (deflection, moment, shear, stress) with `FromResult()` factory
- - [x] `SimulationTimeline` — ordered collection with `FromResult()`, `InterpolateAt()` linear blending
-
-## Phase 4 — ML & Monte Carlo Integration
-
-- [x] `MultiphysicsMonteCarloModel : IMonteCarloModel` — runs solver N times with sampled parameters (*depends on Phase 2; pattern from GIS `PlumeMonteCarloModel`*)
- - [x] `ParameterVariation` — ranges for material properties, BCs, loads
- - [x] Output: scenario matrix, percentile maps
-- [x] `SurrogateTrainer` — trains regression model (Linear/SVR/MLP) on MC scenario data (*depends on MC model*)
- - [x] Predict simulation output without re-running solver
-- [x] `MultiphysicsClusterAnalyzer` — wraps `ClusteringExperiment` on scenario matrix (*parallel with surrogate; pattern from GIS `ScenarioClusterAnalyzer`*)
- - [x] Identifies representative scenarios from MC ensemble
-
-## Phase 5 — Tests & Documentation
-
-- [x] Unit tests with analytical validation
- - [x] Heat plate: steady-state convergence vs known Fourier solutions
- - [x] Pipe flow: velocity profile vs analytical Poiseuille
- - [x] Electric field: capacitor/line charge vs analytical E-field
- - [x] Beam stress: cantilever deflection vs PL³/3EI
- - [x] FEM primitives: stiffness matrix assembly, known beam cases (bar axial, cantilever, simply supported)
- - [x] MC integration: scenario matrix dimensions, clustering produces valid results
- - [x] Export: binary format header parsing, dimension checks
-- [x] Update `Engines/Multiphysics/README.md` — add Multiphysics engine section
-- [x] Update `Numerics/README.md` — add Poisson solver section
-- [x] Update `Physics/README.md` — add EngineeringMaterial section
-
----
-
-## Key Files
-
-### Existing — to modify
-- `Numerics/Numerics/Numerics/FiniteDifference/GridOperators.cs` — add Poisson solver
-- `Numerics/Numerics/Physics/Materials/Materials.cs` — extend factory with Engineering path
-
-### Existing — to reuse as patterns
-- `Numerics/Numerics/Engines/GIS/Simulation/PlumeMonteCarloModel.cs` — MC integration
-- `Numerics/Numerics/Engines/GIS/Analysis/ScenarioClusterAnalyzer.cs` — clustering
-- `Numerics/Numerics/Engines/GIS/Export/UnityBinaryExporter.cs` — binary export format
-- `Numerics/Numerics/Engines/Common/ISimulationEngine.cs` — engine interface
-- `Numerics/Numerics/Physics/HeatExtensions.cs` — HeatEquationRate()
-- `Numerics/Numerics/Physics/FluidExtensions.cs` — NS residuals, pipe flow
-- `Numerics/Numerics/Physics/ElectroMagneticFieldExtensions.cs` — E-field calculations
-
-
-### New — to create
-- `Numerics/Numerics/Numerics/FiniteElement/` — IElement1D, BarElement, BeamElement, Assembler1D, Mesh1D (~6 files)
-- `Numerics/Numerics/Physics/Materials/Engineering/` — EngineeringMaterial, EngineeringLibrary (~2 files)
-- `Numerics/Numerics/Engines/Multiphysics/` — engine, 4 solvers, config, export, snapshots, MC, ML (~15 files)
-
-
-## Decisions
-
-- **FD over LBM** — existing infrastructure sufficient; LBM complexity unjustified for simple geometries
-- **1D FEM primitives in Numerics** — minimal (Element, Assembly, Mesh), beam solver logic in engine
-- **Export like GIS** — UnityBinaryExporter + JSON for web
-- **Engineering materials** — new struct + library extending existing MaterialDescriptor
-- **Pipe flow scope** — 1D Hagen-Poiseuille primary, 2D lid-driven cavity as stretch goal
-- **Beam analysis** — static only initially; dynamic vibration deferred
-- **RL environment** — deferred to follow-up roadmap
\ No newline at end of file
diff --git a/docs/Roadmap-v4.3.md b/docs/Roadmap-v4.3.md
index 27b8697..feedaf3 100644
--- a/docs/Roadmap-v4.3.md
+++ b/docs/Roadmap-v4.3.md
@@ -182,36 +182,75 @@ Punkter som stod i v4.1-scopet och ännu inte är gjorda:
## Implementationsplan — Faser
-### Phase 1 — Rotfinnare
-- [ ] Skapa `Numerics/RootFinding/`-struktur + `RootResult`
-- [ ] Implementera `Bisection`, `Secant`
-- [ ] Implementera `Brent`
-- [ ] Implementera `Newton` med tolerans, maxiter, f′-skydd och analytisk-derivata-overload
-- [ ] `NewtonRaphson` blir wrapper över `Newton` (signatur oförändrad)
-- [ ] Migrera `KeplerOrbit` och `LagrangePoints` till fasaden
-- [ ] Enhetstester: patologiska funktioner, platta derivator, ingen teckenväxling, konvergensrapportering
+### Phase 1 — Rotfinnare ✔ klar
+- [x] Skapa `Numerics/RootFinding/`-struktur + `RootResult`
+- [x] Implementera `Bisection`, `Secant`
+- [x] Implementera `Brent`
+- [x] Implementera `Newton` med tolerans, maxiter, f′-skydd och analytisk-derivata-overload
+- [x] `NewtonRaphson` blir wrapper över `Newton` (signatur oförändrad)
+- [x] Migrera `KeplerOrbit` och `LagrangePoints` till fasaden
+- [x] Enhetstester: patologiska funktioner, platta derivator, ingen teckenväxling, konvergensrapportering
+
+> **Noterat under Phase 1:** `TimeserieValidationTests` har tre fel som *inte* rör rotfinnarna.
+> `TestData/CsvTestDataGenerator` skriver decimaltal med aktuell kultur (`$"{v:F2}"`), så på en svensk
+> maskin blir `6.44` till `6,44` och kolliderar med CSV-avgränsaren. Testerna passerar med
+> `DOTNET_SYSTEM_GLOBALIZATION_INVARIANT=1` och på CI (Linux). Fixen är `CultureInfo.InvariantCulture`
+> i generatorn — ligger utanför v4.3-scopet men bör tas någon gång.
### Phase 2 — Benchmark-baseline
-- [ ] Skapa `Numerics.Benchmarks`-projekt (BenchmarkDotNet), lägg till i `.sln`, exkludera från paketering
-- [ ] Benchmarks för matmul, SpMV, LU/QR/Cholesky, `Solve`
-- [ ] Benchmark för en MLP-träningsepok
-- [ ] Kör och checka in baseline **före** Del 1-migreringen
+- [x] Skapa `Numerics.Benchmarks`-projekt (BenchmarkDotNet), lägg till i `.sln`, exkludera från paketering
+- [x] Benchmarks för matmul, SpMV, LU/QR/Cholesky, `Solve`
+- [x] Benchmark för en MLP-träningsepok
+- [x] Kör och checka in baseline **före** Del 1-migreringen → [baseline-2026-10-01](benchmarks/baseline-2026-10-01.md)
+
+> **Noterat under Phase 2:** Windows app control-policyn (0x800711C7) blockerar den DLL-kopia som
+> BenchmarkDotNet lägger i sin genererade per-benchmark-mapp, så standardtoolchainen (en process per
+> benchmark) ger `No Workload Results` på den här maskinen. Körningarna görs därför med `--inProcess`.
+> Mätningarna är giltiga men något mindre isolerade än med processeparation — jämför alltid baslinje
+> mot omkörning gjord på samma sätt.
### Phase 3 — Migrering till dekompositionerna
-- [ ] Regressionstester som låser nuvarande resultat för de åtta anropsställena
-- [ ] Migrera `MultivariateInterpolation`, `PanelMethod`, `Assembler1D`, `CubicSpline`-fallback, `InferentialStatisticsExtensions`, `DifferentialEquationExtensions` till LU
-- [ ] Migrera `FittingSolver` till QR + verifiera standardfelen mot nuvarande värden
-- [ ] Migrera `KalmanFilter`/`ExtendedKalmanFilter`/`KalmanSmoother` till Cholesky-lösning
-- [ ] Migrera `CoupledOscillators` till `EigenDecomposition`
-- [ ] Migrera `PCA` till `EigenDecomposition`
-- [ ] Låt `EigenValues`/`DominantEigenVector`/`EigenVector` delegera till `EigenDecomposition`
-- [ ] Cacha faktoriseringar där flera högerled löses mot samma matris
-
-### Phase 4 — Städning
-- [ ] `NaiveBayes.NumClasses` sätts i `Fit`
-- [ ] Ta bort tempfilsreferensen i csproj
-- [ ] Ta bort dubblerade roadmaps ur `docs/` respektive `docs/completed/`
-- [ ] Ta bort `SARSA._pendingNextAction`
+- [x] Regressionstester för anropsställena — nya sviter för `PanelMethod` och `PCA`, som saknade
+ täckning helt; övriga sites täcks av befintliga tester
+- [x] Migrera `MultivariateInterpolation`, `Assembler1D`, `CubicSpline`-fallback,
+ `InferentialStatisticsExtensions`, `DifferentialEquationExtensions` till LU
+- [x] Migrera `FittingSolver` till QR + verifiera standardfelen mot nuvarande värden — gjort för
+ `LeastSquaresFitter`, `WeightedLeastSquaresFitter`, `RobustFitter`, `ParameterEstimation` och
+ `NonlinearLeastSquaresFitter`s kovarians. LM-iterationen lämnas på dämpade normalekvationer
+ (dämpningen läggs på `JtJ`s diagonal; QR kräver den utvidgade formuleringen [J; √(λD)], en
+ algoritmändring)
+- [x] Migrera `KalmanFilter`/`ExtendedKalmanFilter`/`KalmanSmoother` till Cholesky-lösning
+- [x] Migrera `CoupledOscillators` till `EigenDecomposition`
+- [x] Migrera `PCA` till `EigenDecomposition`
+- [x] Cacha faktoriseringar där flera högerled löses mot samma matris — gäller Kalman-vinsten,
+ där `Cholesky.Solve(Matrix)` nu löser alla kolumner mot en faktorisering
+
+**Två punkter ströks efter att koden lästs:**
+
+- **`PanelMethod` migrerades inte.** Den har varken tester eller anropare. Karaktäriseringstesterna
+ avslöjade att lösaren ger Cp ≈ 1 på varje panel för en cylinder vid noll anfallsvinkel, där
+ potentialflöde ger ett intervall över [−3, 1]. Om det är en bugg eller om en cylinder ligger
+ utanför dess domän (den påtvingar ett Kutta-villkor vid en bakkant cylindern inte har) är
+ ouppklarat — testet ligger `[Ignore]`-markerat med frågan nedskriven. Utan tester och utan
+ anropare är dess duplicering den minst skadliga, och att utreda lösaren är separat arbete.
+- **`EigenValues`/`DominantEigenVector`/`EigenVector` delegerar inte.** Planpunkten antog att de var
+ duplicerade generella egenlösare. De är i stället lågnoggranna hjälpfunktioner vars publika
+ kontrakt är avrundade heltalskvoter — `Math.Abs(Math.Round(c / min))` — asserterat av befintliga
+ tester (`result[0] == 2`) och konsumerat av `OdeSolver`. Att delegera dem byter ut kontraktet mot
+ normaliserade egenvektorer, vilket är en API-ändring med omvalidering av `OdeSolver`, inte en
+ refaktorering bakom befintligt API. Kräver ett eget beslut.
+
+### Phase 4 — Städning ✔ klar
+- [x] `NaiveBayes.NumClasses` sätts i `Fit` — plus en testsvit, modellen hade ingen
+- [x] Ta bort tempfilsreferensen i csproj
+- [x] Ta bort dubblerade roadmaps ur `docs/` respektive `docs/completed/`
+- [x] Ta bort `SARSA._pendingNextAction`
+
+> **Noterat:** `docs/`-kopian av `AdvancedGameEngineRoadMap` hade blivit helt obockad medan
+> `completed/`-kopian var bockad — den senare stämmer med koden, så dubbletten i `docs/` ströks.
+> `Multiphysics`-paret var omvänt: `docs/`-kopian var den aktuella och fick ersätta den i
+> `completed/`. Två CS0219-varningar återstår i testprojektet (oanvända lokala variabler i
+> `FiniteElementTests` och `SignalProcessingFilterTests`) — de stod inte i scopet och lämnades.
### Phase 5 — Verifiering & release
- [ ] Kör om benchmarks och jämför mot baseline från Phase 2
diff --git a/docs/TerrainSpreadRoadMap.md b/docs/TerrainSpreadRoadMap.md
deleted file mode 100644
index db0238e..0000000
--- a/docs/TerrainSpreadRoadMap.md
+++ /dev/null
@@ -1,319 +0,0 @@
-# Terrain Spread Framework — Wildfire (Rothermel) MVP Roadmap
-
-## Overview
-
-A **terrain-aware spread simulation framework** inside `Engines/GIS/`, starting with **wildfire** as the first scenario. The fire physics are based on the **Rothermel (1972) surface fire spread model** — the standard used by FARSITE, FlamMap, and BehavePlus.
-
-The framework reuses the existing GIS engine infrastructure: `GeoGrid` for the spatial domain, `GridSnapshot` + named layers for per-cell state (fuel, moisture, flame length, burn time), the Monte Carlo / clustering pipeline for stochastic weather ensembles, `ExposurePolygonGenerator` for fire-perimeter extraction, and the export pipeline (GeoJSON, Cesium, Unity) for visualization.
-
-### Architecture
-
-```
-Physics/
-├── Materials/
-│ ├── Chemical/ ← existing (ChemicalSubstance, ChemicalLibrary)
-│ ├── Nuclear/ ← existing (Isotope, IsotopeLibrary, Decay, …)
-│ ├── Engineering/ ← existing (EngineeringMaterial, EngineeringLibrary)
-│ └── Fire/ ← NEW — fuel model data (same pattern as Chemical/Nuclear)
-│ ├── FuelModel.cs Rothermel fuel parameters (Anderson 13)
-│ ├── FuelLibrary.cs Static registry of standard fuel models
-│ └── Enums/
-│ └── FuelModelType.cs Enum: ShortGrass, TimberLitter, …
-├── Environmental/
-│ └── Fire/ ← NEW — fire spread physics (pure math, no grid)
-│ └── RothermelModel.cs Core Rothermel equations (R, IR, φw, φs)
-
-
-Engines/GIS/
-├── Grid/ ← existing (GeoGrid, GridSnapshot, GeoCell)
-├── Terrain/ ← NEW — terrain model + fuel map (grid-aware)
-│ ├── TerrainGrid.cs Elevation surface, slope/aspect
-│ └── FuelMap.cs Per-cell fuel assignment on grid
-├── Spread/ ← NEW — generic spread engine
-│ ├── ISpreadSimulator.cs Interface for any terrain spread model
-│ ├── SpreadResult.cs Timeline of SpreadSnapshot (layers per step)
-│ ├── SpreadSnapshot.cs Per-step cell state (burning, burned, unburned, flame length, ROS)
-│ └── Wildfire/
-│ ├── WildfireSimulator.cs Cell-automaton spread using Physics.Environmental.Fire.RothermelModel
-│ ├── WildfireParameters.cs Runtime config (ignition, wind, moisture)
-│ └── Enums/
-│ └── CellBurnState.cs Unburned, Burning, Burned, Firebreak
-├── Scenario/ ← extend existing fluent API
-│ ├── RiskScenario.cs + ForWildfire() entry point
-│ ├── WildfireScenarioBuilder.cs Fluent configuration for fire scenarios
-│ └── WildfireScenarioResult.cs Fire-specific result (perimeters, area burned, …)
-├── Simulation/ ← existing plume (unchanged)
-├── Analysis/ ← reuse ExposurePolygonGenerator for fire perimeters
-├── Export/ ← extend GeoJSON/Cesium with fire-specific features
-└── RL/ ← future: fire suppression RL environment
-```
-
-### Cross-section Dependencies
-
-```
-Physics.Environmental.Fire.RothermelModel → Physics.Materials.Fire.FuelModel (pure math, no grid deps)
-Physics.Materials.Fire.FuelLibrary → Physics.Materials.Fire.FuelModel
-
-Engines.GIS.WildfireSimulator → Physics.Environmental.Fire.RothermelModel (fire spread physics)
- → Physics.Materials.Fire.FuelModel (fuel data)
- → Engines.GIS.TerrainGrid (slope/aspect)
- → Engines.GIS.FuelMap (fuel per cell, wraps FuelModel)
- → Engines.GIS.GeoGrid (spatial domain, existing)
- → Engines.GIS.GridSnapshot (layers, existing)
-```
-
-This follows the same separation as the plume pipeline:
-- **Materials:** `Physics.Materials.Fire.FuelModel` sits alongside `Physics.Materials.Chemical` and `Physics.Materials.Nuclear`
-- **Physics:** `Physics.Environmental.Fire.RothermelModel` extends the environmental physics alongside `EnvironmentalExtensions.GaussianPlume()`
-- **Engine:** `Engines.GIS.WildfireSimulator` wires everything to the grid, just like `Engines.GIS.PlumeSimulator`
-
-### Rothermel Model Summary
-
-Rate of spread (m/min):
-
-$$R = \frac{I_R \cdot \xi \cdot (1 + \phi_w + \phi_s)}{\rho_b \cdot \varepsilon \cdot Q_{ig}}$$
-
-| Symbol | Meaning |
-|--------|---------|
-| $I_R$ | Reaction intensity (kJ/m²·min) |
-| $\xi$ | Propagating flux ratio |
-| $\phi_w$ | Wind correction factor |
-| $\phi_s$ | Slope correction factor |
-| $\rho_b$ | Ovendry bulk density (kg/m³) |
-| $\varepsilon$ | Effective heating number |
-| $Q_{ig}$ | Heat of pre-ignition (kJ/kg) |
-
-Fuel parameters per model: surface-area-to-volume ratio σ, fuel bed depth δ, ovendry fuel load $w_0$, dead fuel moisture of extinction $M_x$, low heat content $h$.
-
----
-
-## Phase 1 — Fire Physics & Fuel Library (Physics section)
-
-Pure math and data — no grid or simulation. Two namespaces following existing conventions:
-- **Materials:** `CSharpNumerics.Physics.Materials.Fire` (alongside Chemical, Nuclear, Engineering)
-- **Physics:** `CSharpNumerics.Physics.Environmental.Fire` (alongside GaussianPlume)
-
-### Fuel models (Anderson 13) — `Physics/Materials/Fire/`
-
-- [x] `FuelModel` immutable record: `FuelModelType Type`, `string Name`, `double SurfaceAreaToVolumeRatio` (1/m), `double FuelBedDepth` (m), `double OvendryFuelLoad` (kg/m²), `double MoistureOfExtinction` (fraction), `double LowHeatContent` (kJ/kg), `double ParticleDensity` (kg/m³ — default 513 for wood)
-- [x] `FuelModelType` enum for Anderson 13: `ShortGrass (1)`, `TimberGrassUnderstory (2)`, `TallGrass (3)`, `Chaparral (4)`, `Brush (5)`, `DormantBrush (6)`, `SouthernRough (7)`, `ClosedTimberLitter (8)`, `HardwoodLitter (9)`, `TimberLitterUnderstory (10)`, `LightLoggingSlash (11)`, `MediumLoggingSlash (12)`, `HeavyLoggingSlash (13)`, `NoFuel (0)`
-- [x] `FuelLibrary` static class: `Get(FuelModelType)`, `TryGet(type, out fuel)`, `All`, `Register(FuelModel)` — pre-loaded with all 13 Anderson models
-
-### Rothermel core equations — `Physics/Environmental/Fire/RothermelModel.cs`
-
-- [x] `RothermelModel` static class with:
- - `RateOfSpread(FuelModel fuel, double moistureContent, double windSpeed, double slopeRadians)` → `double` (m/min)
- - `ReactionIntensity(FuelModel fuel, double moistureContent)` → `double` IR (kJ/m²·min)
- - `WindFactor(FuelModel fuel, double midflameWindSpeed)` → `double` φw
- - `SlopeFactor(double packingRatio, double slopeRadians)` → `double` φs
- - `PropagatingFluxRatio(FuelModel fuel)` → `double` ξ
- - `HeatOfPreignition(double moistureContent)` → `double` Qig (kJ/kg)
- - `EffectiveHeatingNumber(FuelModel fuel)` → `double` ε
- - `PackingRatio(FuelModel fuel)` → `double` β = ρb/ρp
- - `OptimalPackingRatio(FuelModel fuel)` → `double` β_op
- - `FlameLength(double reactionIntensity, double rateOfSpread)` → `double` (m) — Byram's fireline intensity → flame length correlation
-
-### Tests — Phase 1
-
-- [x] `FuelLibrary` returns all 13 Anderson models, verify Short Grass (σ=3500 1/ft ≈ 11483 1/m, δ=1 ft ≈ 0.305 m)
-- [x] Short Grass (model 1), moisture 0.05, wind 5 mph, flat → R ≈ 23–25 m/min (BehavePlus reference)
-- [x] Chaparral (model 4), moisture 0.10, wind 10 mph, flat → verify against BehavePlus
-- [x] Slope factor: flat terrain → φs = 0
-- [x] Slope factor: 30° slope → φs > 0, increasing with slope
-- [x] Wind factor: zero wind → φw = 0
-- [x] Moisture at extinction → R = 0 (no spread)
-- [x] Moisture above extinction → R = 0 (clamp)
-- [x] Flame length proportional to fireline intensity
-
----
-
-## Phase 2 — Terrain Model & Fuel Map (GIS engine)
-
-Build the elevation surface and per-cell fuel assignment. These are grid-aware wrappers that live in `Engines/GIS/Terrain/` and consume `Physics.Fire.FuelModel`.
-
-### Terrain grid — `Engines/GIS/Terrain/TerrainGrid.cs`
-
-- [x] `TerrainGrid` class wrapping a `GeoGrid` (ground plane, Nz=1) with a `double[] Elevation` array (one height per (ix,iy) cell)
-- [x] `FromFunction(grid, Func elevationFn)` — procedural elevation from f(x,y)
-- [x] `FromArray(grid, double[,] elevation)` — load from 2D array (row = iy, col = ix)
-- [x] `Slope(ix, iy)` — terrain slope in radians using central-difference gradient: $\tan(\theta) = \sqrt{(\partial z/\partial x)^2 + (\partial z/\partial y)^2}$
-- [x] `Aspect(ix, iy)` — downslope direction in radians (0=N, π/2=E, π=S, 3π/2=W)
-- [x] `SlopeInDirection(ix, iy, Vector direction)` — slope component along a given heading (needed for directional Rothermel $\phi_s$)
-
-### Fuel map — `Engines/GIS/Terrain/FuelMap.cs`
-
-- [x] `FuelMap` class: assigns a `Physics.Materials.Fire.FuelModel` per (ix,iy) cell on a `GeoGrid`
- - `SetFuel(ix, iy, FuelModelType)`
- - `SetUniformFuel(FuelModelType)` — fill entire grid
- - `SetFuelByElevation(terrain, ranges)` — assign fuel models to elevation bands
- - `GetFuel(ix, iy)` → `FuelModel`
- - `GetMoisture(ix, iy)` → `double` — per-cell dead fuel moisture content (fraction)
- - `SetMoisture(ix, iy, double)` / `SetUniformMoisture(double)`
-
-### Tests — Phase 2
-
-- [x] `TerrainGrid` slope/aspect on flat surface → slope ≈ 0
-- [x] `TerrainGrid` slope/aspect on known tilted plane → verify against analytical result
-- [x] `SlopeInDirection` matches full slope when direction = aspect, zero when perpendicular
-- [x] `FuelMap` set/get round-trip, uniform fill, moisture defaults
-
----
-
-## Phase 3 — Cell-Automaton Fire Spread Simulator
-
-Wire `Physics.Environmental.Fire.RothermelModel` to the `GeoGrid` via a cellular automaton that propagates fire across the terrain surface.
-
-### Spread engine
-
-- [x] `CellBurnState` enum: `Unburned`, `Burning`, `Burned`, `Firebreak`
-- [x] `ISpreadSimulator` interface:
- ```csharp
- IReadOnlyList Run(GeoGrid grid, TerrainGrid terrain, FuelMap fuelMap,
- TimeFrame timeFrame);
- ```
-- [x] `SpreadSnapshot` class — extends or wraps `GridSnapshot` with:
- - Layer `"burnState"` (0=Unburned, 1=Burning, 2=Burned, 3=Firebreak)
- - Layer `"flameLength"` (metres)
- - Layer `"rateOfSpread"` (m/min)
- - Layer `"burnTime"` (seconds since ignition, 0 if unburned)
- - `BurningCellCount`, `BurnedCellCount`, `BurnedAreaHectares`
-- [x] `WildfireParameters`:
- - `IgnitionPoints` — list of (ix,iy) or world positions
- - `MidflameWindSpeed` (m/s)
- - `WindDirection` (Vector — same convention as PlumeSimulator)
- - `BurnDuration` (seconds — how long a cell stays in Burning state before transitioning to Burned)
- - `SpotFireEnabled` (bool) — future, default false
-- [x] `WildfireSimulator : ISpreadSimulator`
- - 8-neighbour spread (N/S/E/W + diagonals, diagonal distance = step√2)
- - For each Burning cell, compute ROS towards each unburned neighbour:
- 1. Direction from burning cell to neighbour
- 2. Slope in that direction from `TerrainGrid.SlopeInDirection()`
- 3. Wind component along that direction
- 4. Neighbour fuel + moisture → `RothermelModel.RateOfSpread()`
- 5. Travel time = distance / ROS
- 6. If accumulated time ≥ travel time → ignite neighbour
- - Time step loop: advance time, update burning→burned transitions, ignite reachable neighbours
- - Produce one `SpreadSnapshot` per time step
-
-### Tests — Phase 3
-
-- [x] Single ignition on flat uniform Short Grass with no wind → near-circular spread pattern
-- [x] Ignition with steady wind → elliptical spread (elongated downwind)
-- [x] Uphill slope accelerates spread, downhill decelerates
-- [x] `NoFuel` cells block fire (act as firebreaks)
-- [x] `Firebreak` cell state blocks spread
-- [x] `BurnedAreaHectares` increases monotonically
-- [x] Zero ROS at high moisture (saturated fuel) → fire does not spread
-
----
-
-## Phase 4 — Fluent API & Scenario Integration
-
-Expose wildfire through the same fluent builder pattern as the plume scenario, including Monte Carlo for weather uncertainty.
-
-### Builder
-
-- [x] `RiskScenario.ForWildfire()` → returns `WildfireScenarioBuilder`
-- [x] `WildfireScenarioBuilder` fluent chain:
- ```csharp
- RiskScenario
- .ForWildfire()
- .WithTerrain(terrainGrid)
- .WithFuel(fuelMap)
- .WithIgnition(ix, iy) // or .WithIgnition(position)
- .WithWind(speed, direction)
- .WithMoisture(0.08) // global dead fuel moisture
- .OverGrid(grid)
- .OverTime(0, 7200, 60) // 2 hours, 1-min steps
- .RunSingle(); // → WildfireScenarioResult
- ```
-- [x] `WildfireScenarioResult` — like `ScenarioResult` but fire-specific:
- - `Snapshots` — `List`
- - `FinalBurnedArea` (hectares)
- - `MaxFlameLength` (metres)
- - `FirePerimeters` — list of `ExposurePolygon` (one per time step, from burn boundary)
- - `GenerateFirePerimeter(timeIndex)` → `ExposurePolygon` via `ExposurePolygonGenerator` on `burnState ≥ 1`
-
-### Monte Carlo
-
-- [x] `WildfireVariation` — stochastic parameter ranges:
- - Wind speed range
- - Wind direction jitter
- - Moisture content range
- - Ignition location offset radius
-- [x] `.WithVariation(v => v.WindSpeed(3, 8).Moisture(0.04, 0.12))`
-- [x] `.RunMonteCarlo(iterations)` → `WildfireMonteCarloResult`
- - Per-cell burn probability across all iterations
- - Mean / max burned area statistics
-- [x] `.AnalyzeWith(clustering)` — reuse existing `ScenarioClusterAnalyzer` on the burn-probability matrix
-- [x] `.Build()` → `WildfireScenarioResult` with probability-weighted outputs
-
-### Tests — Phase 4
-
-- [x] Fluent API deterministic round-trip: build → run → inspect burned area
-- [x] Monte Carlo 20 iterations: burn probability ∈ [0, 1] for all cells
-- [x] Clustering identifies distinct fire spread regimes (high-wind vs. low-wind clusters)
-- [x] `GenerateFirePerimeter()` returns valid polygon enclosing burned cells
-
----
-
-## Phase 5 — Export & Visualization
-
-Extend the existing export pipeline with fire-specific outputs.
-
-### GeoJSON
-
-- [x] Point features with `burnState`, `flameLength`, `rateOfSpread` properties
-- [x] Fire perimeter as `Polygon` geometry per time step (via `ExposurePolygonGenerator`)
-- [x] Burn probability heatmap export (from Monte Carlo)
-
-### Cesium (CZML)
-
-- [x] Time-dynamic fire perimeter polygons (animated spread)
-- [x] Colour ramp: red (burning) → grey (burned) → green (unburned)
-
-### Unity binary
-
-- [x] Extend `UnityBinaryExporter` to write fire layers (burnState, flameLength)
-
-### Tests — Phase 5
-
-- [x] GeoJSON export produces valid FeatureCollection with fire perimeter polygons
-- [x] CZML export has time intervals matching simulation time steps
-- [x] Binary round-trip: write → read → verify layer values
-
----
-
-## Phase 6 — Documentation & Polish
-
-- [x] Update `Physics/README.md` with Fire section:
- - RothermelModel overview and equations
- - FuelModel / FuelLibrary / FuelModelType reference (Physics.Materials.Fire)
- - Standalone usage examples (ROS calculation without grid)
-- [x] Update `Engines/GIS/README.md` with Wildfire / Terrain Spread section:
- - TerrainGrid, FuelMap overview
- - WildfireSimulator usage (consuming Physics.Environmental.Fire + Physics.Materials.Fire)
- - Fluent API examples
- - Monte Carlo fire probability workflow
- - Export examples
-- [x] Add code samples for common scenarios:
- - Flat grassland fire
- - Mountainous terrain with mixed fuel
- - MC ensemble with wind uncertainty
-- [x] Validate all public types have XML doc summaries
-- [x] Final test pass — all wildfire tests green
-
----
-
-## Future Extensions (Out of Scope for MVP)
-
-These are **not** part of the MVP but inform the architecture to keep extensibility open:
-
-- **Spot fires** — ember transport via lofting model (wind + convection column)
-- **Crown fire** — Van Wagner (1977) crown fire initiation + spread
-- **Scott & Burgan 40 fuel models** — expanded fuel library
-- **Weather timeline** — time-varying wind speed/direction/moisture during simulation
-- **Suppression modelling** — firebreaks, water drops, crew lines (RL environment candidate)
-- **Additional spread scenarios** — flood/inundation, landslide, disease/pest, pollutant propagation (all reuse `ISpreadSimulator` + `TerrainGrid`)
-- **DEM import** — read GeoTIFF / ASCII grid elevation data
-- **Fuel map import** — read LANDFIRE / Corine raster land cover → fuel model mapping
diff --git a/docs/benchmarks/baseline-2026-10-01.md b/docs/benchmarks/baseline-2026-10-01.md
new file mode 100644
index 0000000..31a6539
--- /dev/null
+++ b/docs/benchmarks/baseline-2026-10-01.md
@@ -0,0 +1,114 @@
+# Prestandabaslinje — 2026-10-01
+
+Baslinjen som v4.3:s migrering och v4.4:s prestandaarbete mäts emot. Tagen **före** Del 1
+(migreringen av anropsställena till dekompositionerna), enligt principen *mät före optimering*.
+
+| | |
+|---|---|
+| **Commit** | `7832fe8` (gren `feat/root-finding`) |
+| **Maskin** | Intel Core i5-14400 2.50GHz, 16 logiska / 10 fysiska kärnor, Windows 11 25H2 |
+| **Runtime** | .NET 10.0.9, X64 RyuJIT x86-64-v3, Release |
+| **Verktyg** | BenchmarkDotNet 0.15.8, `InProcessEmitToolchain` |
+| **Körtid** | 7 min 56 s, 22 mätpunkter |
+
+> **Om toolchainen:** Windows app control-policyn (0x800711C7) blockerar den DLL-kopia som
+> BenchmarkDotNet lägger i sin genererade per-benchmark-mapp, så standardtoolchainen ger
+> `No Workload Results` på den här maskinen. Körningen gjordes därför med `--inProcess`.
+> Siffrorna är giltiga men mindre isolerade än med processeparation — **omkörningar måste göras
+> på samma sätt för att vara jämförbara.**
+
+Reproducera:
+
+```powershell
+dotnet run -c Release --project Numerics/Numerics.Benchmarks -- --filter * --inProcess
+```
+
+---
+
+## Täta kärnor
+
+| Metod | N | Tid | Allokerat |
+|-------|---|----:|----------:|
+| Matrix × Matrix | 64 | 425.9 µs | 64.1 KB |
+| Matrix × Vector | 64 | 2.9 µs | 1.1 KB |
+| LU-faktorisering | 64 | 106.4 µs | 32.4 KB |
+| QR-faktorisering | 64 | 218.4 µs | 32.6 KB |
+| Cholesky-faktorisering | 64 | 49.9 µs | 32.1 KB |
+| LU faktorisera + lös | 64 | 112.2 µs | 33.4 KB |
+| LU lös (återanvänd faktorisering) | 64 | **6.6 µs** | 1.1 KB |
+| `Matrix.Inverse` | 64 | 538.6 µs | 194.0 KB |
+| Matrix × Matrix | 256 | 29,737.6 µs | 1,024.3 KB |
+| Matrix × Vector | 256 | 40.5 µs | 4.1 KB |
+| LU-faktorisering | 256 | 6,193.4 µs | 514.3 KB |
+| QR-faktorisering | 256 | 14,232.6 µs | 514.9 KB |
+| Cholesky-faktorisering | 256 | 2,778.2 µs | 512.5 KB |
+| LU faktorisera + lös | 256 | 6,782.0 µs | 518.3 KB |
+| LU lös (återanvänd faktorisering) | 256 | **118.6 µs** | 4.1 KB |
+| `Matrix.Inverse` | 256 | 35,459.3 µs | 3,081.3 KB |
+
+## Glesa kärnor
+
+5-punkts-Laplacian över inre noder på ett `Grid`×`Grid`-nät.
+
+| Metod | Grid | Obekanta | Tid | Allokerat |
+|-------|------|---------:|----:|----------:|
+| SpMV | 32 | 900 | 3.8 µs | 14.1 KB |
+| PCG-lösning | 32 | 900 | 398.5 µs | 1,608.5 KB |
+| SpMV | 128 | 15,876 | 71.8 µs | 248.4 KB |
+| PCG-lösning | 128 | 15,876 | 33,588.1 µs | 117,750.9 KB |
+
+## Maskininlärning
+
+| Metod | Samples | Tid | Allokerat |
+|-------|---------|----:|----------:|
+| MLP-regressor, en träningsepok | 256 | 2.12 ms | 11.13 MB |
+| MLP-regressor, en träningsepok | 2048 | 17.39 ms | 88.75 MB |
+
+---
+
+## Observationer
+
+**1. Att återanvända en faktorisering är 57× snabbare än att faktorisera om.** Vid N=256 kostar
+`LU faktorisera + lös` 6,782 µs medan enbart `Solve` mot en färdig faktorisering kostar 118.6 µs.
+Det här är hela argumentet för v4.3:s regel *faktorisera en gång*: överallt där flera högerled
+löses mot samma matris — `PanelMethod` över anfallsvinklar, FEM-assemblern, Kalman-looparna —
+ligger en storleksordning och väntar.
+
+**2. `Matrix.Inverse` är den dyraste operation som mätts.** 35.5 ms vid N=256, alltså 5× kostnaden
+för en komplett LU-faktorisering med lösning, och 6× allokeringen. Varje ställe som skriver en
+lösning som `A.Inverse() * b` betalar det priset i onödan. Det gäller särskilt
+`KalmanFilter`/`ExtendedKalmanFilter`/`KalmanSmoother`, som alla gör `S.Inverse()` — se
+begränsning 4 i [Roadmap-v4.3](../Roadmap-v4.3.md).
+
+**3. Varje `Matrix` allokerar en andra n×n-array som nästan aldrig används.** En 256×256-matmul
+allokerar 1,024 KB, men resultatdatan är bara 512 KB. Resten är `identity`-fältet, som
+`Matrix`-konstruktorn fyller i vid varje konstruktion oavsett om någon läser det. Det dubblar
+minneskostnaden för varenda matris i biblioteket. Kandidat för v4.4 — men notera att `Matrix` är
+en `struct` med publika fält, så en ändring måste göras varsamt.
+
+**4. `new VectorN(double[])` kopierar sin indata.** SpMV bygger en resultatarray och skickar den
+till konstruktorn, som klonar den: 248 KB allokerat för en vektor vars data är 124 KB. Samma
+dubblering drabbar varje vektorproducerande operation i biblioteket.
+
+**5. Matrismultiplikationen är en naiv trippelloop.** 4× dimensionen ger 70× tiden (425.9 µs →
+29,737.6 µs), mot 64× för ren O(n³) — resten är cachemissar. Det är det primära målet för
+SIMD-arbetet i v4.4, och Gen2-kolumnen visar att 256×256-resultaten går direkt till LOH.
+
+**6. Träningsloopen allokerar 88.75 MB per epok** vid 2048 samples — 44 KB per sample. Det är
+referenspunkten för de allokeringsfria träningslooparna i v4.4.
+
+**7. PCG-lösningen allokerar 117 MB för ett enda anrop** på 15,876 obekanta. Observation 4
+förklarar hälften; resten är att varje CG-iteration allokerar nya vektorer för residual,
+sökriktning och prekonditionerad residual. En `Span`-baserad eller buffertåteranvändande variant
+skulle ta bort nästan allt.
+
+---
+
+## Vad som ska hända med de här siffrorna
+
+Del 1-migreringen i v4.3 **ska inte göra något långsammare**. Den flyttar anropsställen till samma
+dekompositioner som mäts här, så täta kärnor förväntas ligga stilla medan anropsställen som
+tidigare faktoriserade om per lösning blir snabbare. Kör om och jämför innan v4.3 taggas.
+
+Observation 3, 4, 5, 6 och 7 är *inte* v4.3-arbete — de är underlaget för v4.4 och noteras här för
+att de upptäcktes nu.
diff --git a/docs/completed/Multiphysics-Roadmap.md b/docs/completed/Multiphysics-Roadmap.md
index 174154c..b44165c 100644
--- a/docs/completed/Multiphysics-Roadmap.md
+++ b/docs/completed/Multiphysics-Roadmap.md
@@ -25,12 +25,12 @@ Lattice Boltzmann Methods were evaluated and **rejected**. FD (Grid2D, GridOpera
- [x] Add iterative Poisson solver (Gauss-Seidel) to `GridOperators` for ∇²φ = f on Grid2D with Dirichlet BC — `GridOperators.SolvePoisson2D()`
- [x] Add `SolidExtensions` to `Physics/` — Hooke's law, second moment of area, Euler-Bernoulli beam equation, analytical deflections
-- [ ] Add minimal 1D FEM primitives in `Numerics/Numerics/Numerics/FiniteElement/` (*deferred — analytical solvers used for now*)
- - [ ] `IElement1D` interface — shape functions, local stiffness, local load
- - [ ] `BarElement` (2-node, linear) — axial stress/strain
- - [ ] `BeamElement` (2-node, Hermite cubic) — Euler-Bernoulli bending
- - [ ] `Assembler1D` — local→global stiffness matrix assembly into `Matrix`
- - [ ] `Mesh1D` — 1D mesh from interval subdivision (nodes + elements)
+- [x] Add minimal 1D FEM primitives in `Numerics/Numerics/Numerics/FiniteElement/`
+ - [x] `IElement1D` interface — shape functions, local stiffness, local load
+ - [x] `BarElement` (2-node, linear) — axial stress/strain
+ - [x] `BeamElement` (2-node, Hermite cubic) — Euler-Bernoulli bending
+ - [x] `Assembler1D` — local→global stiffness matrix assembly, solved via `LuDecomposition`
+ - [x] `Mesh1D` — 1D mesh from interval subdivision (nodes + elements)
- [x] Add `EngineeringMaterial` to `Physics/Materials/Engineering/`
- [x] `EngineeringMaterial` immutable struct: ThermalConductivity, SpecificHeat, Density, DynamicViscosity, ElectricPermittivity, YoungsModulus, PoissonsRatio
- [x] `EngineeringLibrary` with common materials: steel, aluminum, copper, water, air, concrete, glass
@@ -98,7 +98,7 @@ Lattice Boltzmann Methods were evaluated and **rejected**. FD (Grid2D, GridOpera
- [x] Pipe flow: velocity profile vs analytical Poiseuille
- [x] Electric field: capacitor/line charge vs analytical E-field
- [x] Beam stress: cantilever deflection vs PL³/3EI
- - [ ] FEM primitives: stiffness matrix assembly, known beam cases (*deferred — no FEM implemented*)
+ - [x] FEM primitives: stiffness matrix assembly, known beam cases (bar axial, cantilever, simply supported)
- [x] MC integration: scenario matrix dimensions, clustering produces valid results
- [x] Export: binary format header parsing, dimension checks
- [x] Update `Engines/Multiphysics/README.md` — add Multiphysics engine section