From 2de09b469ea1149eda9ead85163accce9e408276 Mon Sep 17 00:00:00 2001 From: backlundtransform Date: Wed, 30 Sep 2026 14:24:40 +0200 Subject: [PATCH 1/8] feat(numerics): add the root-finding module The library had exactly one root finder, NumericExtensions.NewtonRaphson, and it could not be trusted. It ran a fixed 100 iterations with no convergence test and no early exit, divided by the derivative without checking it for zero - yielding a silent Inf or NaN where the tangent went flat - took no tolerance or iteration budget, and gave the caller no way to tell a converged answer from a diverged one. Each of those 100 iterations called DerivativeExtensions.Derivate, which allocates a Pascal matrix per call, so an easy root cost 100 matrix allocations and 200 function evaluations no matter how early it actually converged. Add Numerics/RootFinding/ with Bisection, Secant, Brent and a real Newton, each returning a RootResult carrying the estimate, whether it converged, the iteration count and the residual, so failure is reported rather than hidden. Brent is the default for a bracketed search and is exposed as the FindRoot extension. Newton stops on a vanishing or non-finite slope and on a non-finite iterate instead of iterating on NaN, and takes an analytic derivative through an overload. NewtonRaphson keeps its signature and delegates to the new Newton. Its default derivative reproduces Derivate at order 1 exactly - the same backward difference over the same 1e-5 step - so existing results are unchanged, but without the per-iteration allocation and with an early exit. Migrate the two hand-rolled Newton loops in the physics module, KeplerOrbit.TrueAnomaly and LagrangePoints.SolveCollinear, onto the shared implementation. Both were already correct; this removes the duplication. Also make NumericsTests.FindRoots assert the root it computes - it asserted Math.Abs(2) == 2 and never looked at the result. Suite: 1570 passing. The three failures in TimeserieValidationTests are unrelated and predate this branch: TestData/CsvTestDataGenerator formats decimals with the current culture, so on a Swedish machine 6.44 is written as 6,44 and collides with the CSV delimiter. They pass under DOTNET_SYSTEM_GLOBALIZATION_INVARIANT=1 and on CI. --- Numerics/NumericTest/NumericsTests.cs | 2 +- Numerics/NumericTest/RootFindingTests.cs | 346 +++++++++++++ .../Numerics/Numerics/NumericExtensions.cs | 25 +- .../Numerics/RootFinding/RootFinder.cs | 460 ++++++++++++++++++ .../RootFinding/RootFindingExtensions.cs | 34 ++ .../Numerics/RootFinding/RootResult.cs | 61 +++ .../Numerics/Physics/Astro/KeplerOrbit.cs | 18 +- .../Physics/Gravitation/LagrangePoints.cs | 34 +- docs/Roadmap-v4.3.md | 22 +- 9 files changed, 953 insertions(+), 49 deletions(-) create mode 100644 Numerics/NumericTest/RootFindingTests.cs create mode 100644 Numerics/Numerics/Numerics/RootFinding/RootFinder.cs create mode 100644 Numerics/Numerics/Numerics/RootFinding/RootFindingExtensions.cs create mode 100644 Numerics/Numerics/Numerics/RootFinding/RootResult.cs 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/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/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/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/docs/Roadmap-v4.3.md b/docs/Roadmap-v4.3.md index 27b8697..d50133f 100644 --- a/docs/Roadmap-v4.3.md +++ b/docs/Roadmap-v4.3.md @@ -182,14 +182,20 @@ 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 From 7a51322ce4d09938668d8d33eff32d2863a1878f Mon Sep 17 00:00:00 2001 From: backlundtransform Date: Thu, 1 Oct 2026 15:25:02 +0200 Subject: [PATCH 2/8] docs(numerics): document the root-finding module in the section README AGENTS.md requires the section README to be updated when public API lands in that section, and the root-finding module added a namespace without one. Add a Root Finding section covering FindRoot, the four methods and when each applies, RootResult and EnsureConverged, and point the existing Newton-Raphson bullet at it. --- Numerics/Numerics/Numerics/README.md | 59 ++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) 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 From 7832fe80377957d0e07999f37127bd949cbc61b1 Mon Sep 17 00:00:00 2001 From: backlundtransform Date: Thu, 1 Oct 2026 15:29:32 +0200 Subject: [PATCH 3/8] phase 2 --- .gitignore | 4 + .../DenseLinearAlgebraBenchmarks.cs | 118 ++++++++++++++++++ .../MachineLearningBenchmarks.cs | 72 +++++++++++ .../Numerics.Benchmarks.csproj | 27 ++++ Numerics/Numerics.Benchmarks/Program.cs | 12 ++ Numerics/Numerics.Benchmarks/README.md | 65 ++++++++++ .../SparseLinearAlgebraBenchmarks.cs | 99 +++++++++++++++ Numerics/Numerics.sln | 14 +++ docs/Roadmap-v4.3.md | 14 ++- 9 files changed, 421 insertions(+), 4 deletions(-) create mode 100644 Numerics/Numerics.Benchmarks/DenseLinearAlgebraBenchmarks.cs create mode 100644 Numerics/Numerics.Benchmarks/MachineLearningBenchmarks.cs create mode 100644 Numerics/Numerics.Benchmarks/Numerics.Benchmarks.csproj create mode 100644 Numerics/Numerics.Benchmarks/Program.cs create mode 100644 Numerics/Numerics.Benchmarks/README.md create mode 100644 Numerics/Numerics.Benchmarks/SparseLinearAlgebraBenchmarks.cs 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/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/docs/Roadmap-v4.3.md b/docs/Roadmap-v4.3.md index d50133f..8e3ee9b 100644 --- a/docs/Roadmap-v4.3.md +++ b/docs/Roadmap-v4.3.md @@ -198,10 +198,16 @@ Punkter som stod i v4.1-scopet och ännu inte är gjorda: > 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 + +> **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 From 3024a82bb0d102727c63b7db00358ad7cacf913e Mon Sep 17 00:00:00 2001 From: backlundtransform Date: Thu, 1 Oct 2026 15:34:29 +0200 Subject: [PATCH 4/8] docs(benchmarks): record the v4.3 performance baseline Measured before the Del 1 migration, so the migration can be shown not to cost anything and the v4.4 performance work has something to beat. The numbers that bear on v4.3 directly: reusing an LU factorization is 57x faster than refactoring per solve at N=256 (119us against 6782us), which is the whole argument for the "factor once" rule; and Matrix.Inverse is the most expensive operation measured at 35.5ms, five times a full LU factor-and-solve, which is what the Kalman filters pay today by forming S.Inverse() explicitly. Four findings belong to v4.4 rather than v4.3 but were discovered here, so they are written down: every Matrix allocates a second n*n array for the identity field that almost nothing reads, doubling the memory cost of every matrix; new VectorN(double[]) clones its input, so SpMV allocates its result twice; the triple-loop multiply grows 70x for a 4x dimension against 64x for pure O(n^3), the rest being cache misses; and one MLP training epoch allocates 88.75 MB over 2048 samples. Runs use --inProcess: the Windows application control policy blocks the DLL copy BenchmarkDotNet places in its generated per-benchmark folder, so the default toolchain yields no results on this machine. Re-runs must match for the comparison to mean anything. --- docs/Roadmap-v4.3.md | 2 +- docs/benchmarks/baseline-2026-10-01.md | 114 +++++++++++++++++++++++++ 2 files changed, 115 insertions(+), 1 deletion(-) create mode 100644 docs/benchmarks/baseline-2026-10-01.md diff --git a/docs/Roadmap-v4.3.md b/docs/Roadmap-v4.3.md index 8e3ee9b..8125ef6 100644 --- a/docs/Roadmap-v4.3.md +++ b/docs/Roadmap-v4.3.md @@ -201,7 +201,7 @@ Punkter som stod i v4.1-scopet och ännu inte är gjorda: - [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 +- [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 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. From 3f81b1276e4f6e5ee3439d8391e974620fc8bf27 Mon Sep 17 00:00:00 2001 From: backlundtransform Date: Fri, 2 Oct 2026 09:42:17 +0200 Subject: [PATCH 5/8] refactor: route the dense solve sites onto LuDecomposition v4.1 built the decompositions but left every caller solving systems its own way. Five of those sites now go through LuDecomposition instead of their own Gaussian elimination: the GaussElimination extension itself, the RBF weight solve in MultivariateInterpolation, the polynomial regression solve in InferentialStatisticsExtensions, the dense not-a-knot fallback in CubicSplineInterpolation, and the global system in Assembler1D. The tridiagonal branch of the spline keeps the Thomas algorithm, which is O(n) against LU's O(n^3) and must not be replaced by it. GaussElimination was not merely duplicated, it was wrong. Its back substitution assigned the pivot A[i,i] to the solution component whenever the numerator came out exactly zero, so any system with a zero in its solution was answered incorrectly and silently: the diagonal system 2x=2, 3y=0, 4z=4 returned [1, 3, 1] instead of [1, 0, 1]. It now returns the right answer, and no longer modifies its arguments in place. The existing tests assert exact equality and still hold. Add LuDecomposition.SmallestPivotMagnitude so callers can detect a near-singular matrix rather than only an exactly singular one, and express IsSingular in terms of it. Without this the RBF solve would have lost its near-singularity guard and reported a bare "not invertible" in place of advice to change epsilon or kernel. Add characterisation tests for PanelMethod, which had neither tests nor callers. They assert invariants rather than physical correctness, because the one check against a closed-form solution does not pass: for a 200-panel cylinder at zero incidence the solver returns Cp ~= 1 on every panel, where potential flow gives a distribution over [-3, 1]. Whether the solver is wrong or a seamed cylinder is outside its domain is unresolved, so that test is left ignored with the question recorded, and PanelMethod itself is not migrated. Suite: 1580 passing, 1 ignored. --- Numerics/NumericTest/PanelMethodTests.cs | 206 ++++++++++++++++++ .../DifferentialEquationExtensions.cs | 72 +----- .../Numerics/FiniteElement/Assembler1D.cs | 45 +--- .../Interpolation/CubicSplineInterpolation.cs | 41 +--- .../MultivariateInterpolation.cs | 61 ++---- .../Decompositions/LuDecomposition.cs | 24 +- .../InferentialStatisticsExtensions.cs | 50 ++--- 7 files changed, 289 insertions(+), 210 deletions(-) create mode 100644 Numerics/NumericTest/PanelMethodTests.cs 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/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/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; } /// From a7eddb80af64ed2fab61cebeff18b11bb38492aa Mon Sep 17 00:00:00 2001 From: backlundtransform Date: Fri, 2 Oct 2026 10:24:11 +0200 Subject: [PATCH 6/8] refactor: solve via Cholesky and EigenDecomposition in the remaining sites Three migrations, continuing the work of routing callers onto the decompositions v4.1 built. The Kalman filters formed an explicit inverse of the innovation covariance to compute the gain. S is symmetric positive definite, so K = P*Ht*S^-1 is now solved as S*Kt = (P*Ht)t instead. The baseline measured Matrix.Inverse at 35.5ms for a 256x256 matrix against 6.8ms for a full LU factor-and-solve, and an explicit inverse can also break the symmetry the filter depends on. Add Matrix.SolveSymmetricPositiveDefinite, which uses Cholesky and falls back to LU if the matrix is not positive definite, since a recursively updated covariance can drift out of definiteness in floating point. It lives in LinearAlgebra rather than in StateEstimation because that is the section that owns the abstraction. All four inverse sites across KalmanFilter, ExtendedKalmanFilter and KalmanSmoother are converted. CoupledOscillators had its own private Jacobi eigensolver, duplicating EigenDecomposition. Its symmetrised dynamical matrix is symmetric by construction, so the shared decomposition takes its symmetric path and provides the same contract: ascending real eigenvalues, orthonormal eigenvectors as columns. PCA found its eigenpairs by power iteration with deflation. It now uses EigenDecomposition for both the covariance and the dual Gram path. This is an accuracy fix, not only de-duplication: two of the sixteen new PCA tests fail against the old implementation, with score columns correlated to 2.1e-5 and reconstruction off by 7.4e-8, which is what deflation error at a 1e-8 tolerance leaves behind. PowerIteration is deleted. MaxIterations, Tolerance and Seed no longer affect anything but are kept so existing calls and grid-search parameter dictionaries keep working; they are documented as inert. PCA had no tests at all before this. The new suite asserts properties of the decomposition rather than one implementation's output. Two planned items are dropped, with reasons recorded in the roadmap: PanelMethod is not migrated, and the legacy EigenValues, EigenVector and DominantEigenVector helpers do not delegate, because their published contract is rounded integer ratios that OdeSolver and existing tests depend on. Changing it is an API decision, not a refactor. Verification: the solution builds clean in Release on all three target frameworks, and the Kalman (6) and CoupledOscillators (45) tests passed before the Windows application control policy began blocking the test assembly continuously. The full suite has NOT been run against these three changes. Run it, or let CI do it, before merging. --- Numerics/NumericTest/PcaTests.cs | 312 ++++++++++++++++++ .../DimensionalityReduction/Algorithms/PCA.cs | 110 ++---- .../MatrixDecompositionExtensions.cs | 19 ++ .../Oscillations/CoupledOscillators.cs | 119 +------ .../StateEstimation/ExtendedKalmanFilter.cs | 6 +- .../StateEstimation/KalmanFilter.cs | 8 +- .../StateEstimation/KalmanSmoother.cs | 12 +- docs/Roadmap-v4.3.md | 31 +- 8 files changed, 421 insertions(+), 196 deletions(-) create mode 100644 Numerics/NumericTest/PcaTests.cs 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/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/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/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/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/Roadmap-v4.3.md b/docs/Roadmap-v4.3.md index 8125ef6..e432b15 100644 --- a/docs/Roadmap-v4.3.md +++ b/docs/Roadmap-v4.3.md @@ -210,14 +210,31 @@ Punkter som stod i v4.1-scopet och ännu inte är gjorda: > 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 +- [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 - [ ] 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 +- [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 - [ ] `NaiveBayes.NumClasses` sätts i `Fit` From 2b35f4b9ef9d0c4295ab31948352b728604c1194 Mon Sep 17 00:00:00 2001 From: backlundtransform Date: Tue, 6 Oct 2026 16:45:27 +0200 Subject: [PATCH 7/8] refactor(statistics): solve the fitters by QR instead of the normal equations Every linear fitter built the design matrix and then immediately reduced it to XtX before solving. That squares the condition number of X, and a Vandermonde design - any polynomial fit - is poorly conditioned to begin with. The solve now happens on the design matrix itself. Measured on a degree-5 fit over [1, 2], by running this branch's new test against the previous implementation: the normal equations recover the coefficients to a worst-case error of 5.6e-5, QR reaches 7.2e-11. That is close to six decimal digits, and it is why the change is worth making rather than merely tidy. Add two primitives to FittingSolver: SolveLeastSquares, and SolveWeightedLeastSquares which scales each row by the square root of its weight so the weighted problem becomes an ordinary one. Both return (XtX)^-1 alongside the coefficients, recovered from the triangular factor as R^-1 R^-T, so the Gram matrix is never formed or inverted for standard errors either. LeastSquaresFitter, WeightedLeastSquaresFitter, RobustFitter (initial OLS, every IRLS step, and its covariance) and ParameterEstimation all move over, as does the covariance estimate in NonlinearLeastSquaresFitter. The Levenberg-Marquardt iteration itself does not move. Marquardt damping is applied to JtJ's diagonal, which has no design-matrix equivalent short of the augmented formulation [J; sqrt(lambda*D)] - a change to the algorithm, not to the solver, and so a separate decision. It is noted in the code, the roadmap and the section README. Rank deficiency is now detected by QrDecomposition.IsFullRank, which uses a relative threshold, rather than by an absolute 1e-14 pivot test. The exception type is unchanged, so NonlinearLeastSquaresFitter's catch still breaks its loop as before. Add FittingConditioningTests, which pin the accuracy down with a tolerance of 1e-9 - comfortably met by QR, far out of reach of the normal equations - so a revert to XtX fails loudly instead of quietly losing digits. Suite: 1602 passing, 1 ignored, verified in Release configuration. --- .../NumericTest/FittingConditioningTests.cs | 195 ++++++++++++++++++ .../Statistics/Fitting/FittingSolver.cs | 105 ++++++++++ .../Statistics/Fitting/LeastSquaresFitter.cs | 8 +- .../Fitting/NonlinearLeastSquaresFitter.cs | 9 +- .../Statistics/Fitting/ParameterEstimation.cs | 6 +- .../Statistics/Fitting/RobustFitter.cs | 19 +- .../Fitting/WeightedLeastSquaresFitter.cs | 8 +- Numerics/Numerics/Statistics/README.md | 12 +- docs/Roadmap-v4.3.md | 6 +- 9 files changed, 340 insertions(+), 28 deletions(-) create mode 100644 Numerics/NumericTest/FittingConditioningTests.cs 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/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/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/docs/Roadmap-v4.3.md b/docs/Roadmap-v4.3.md index e432b15..6e55216 100644 --- a/docs/Roadmap-v4.3.md +++ b/docs/Roadmap-v4.3.md @@ -214,7 +214,11 @@ Punkter som stod i v4.1-scopet och ännu inte är gjorda: täckning helt; övriga sites täcks av befintliga tester - [x] Migrera `MultivariateInterpolation`, `Assembler1D`, `CubicSpline`-fallback, `InferentialStatisticsExtensions`, `DifferentialEquationExtensions` till LU -- [ ] Migrera `FittingSolver` till QR + verifiera standardfelen mot nuvarande värden +- [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` From 5cacb9dc96decc67af7692fd69f84d0e54645b85 Mon Sep 17 00:00:00 2001 From: backlundtransform Date: Thu, 8 Oct 2026 08:27:29 +0200 Subject: [PATCH 8/8] chore: clear the cleanup backlog carried since the v4.1 scope Four items that had been on the list since v4.1 and never done. NaiveBayes.NumClasses threw NotImplementedException - the only such member left in the library, and the only IClassificationModel that did not provide it. It is now set during Fit as the highest label plus one, matching DecisionTree, RandomForest, KernelSVC, KNearestNeighbors and LinearSVC, so callers can use it to size label-indexed arrays. The model had no tests at all, so add a suite covering that, separation of two Gaussian blobs, classification of unseen points, the unfitted-model guard, and the hyperparameters-only Clone convention. Remove the csproj reference excluding NumericExtensions.cs~RF43b9e4.TMP. The file has not existed for a long time. Remove the unused SARSA._pendingNextAction field, which produced a CS0169 warning on every target framework. Library warnings drop from three to zero; the two that remain are unused locals in the test project and were not in scope. Consolidate four roadmaps that existed in both docs/ and docs/completed/. Three resolved straightforwardly, but the pairs were not simply duplicates and the choice mattered: the docs/ copy of AdvancedGameEngineRoadMap had become entirely unchecked while the completed/ copy was checked, and the completed/ copy is the one that matches the code - the Aerodynamics models are present, and AircraftState and FlightDynamicsEngine are absent only because Engines.Game was extracted into its own package. Multiphysics was the reverse: the docs/ copy recorded the 1-D FEM primitives as done, which they are, while the completed/ copy still called them deferred. That copy is replaced by the accurate one, with its Assembler1D line corrected to say LuDecomposition rather than Gaussian elimination, which this branch changed. Suite: 1608 passing, 1 ignored, verified in Release configuration. --- Numerics/NumericTest/NaiveBayesTests.cs | 111 ++++++ Numerics/Numerics/CSharpNumerics.csproj | 6 +- .../ML/Models/Classification/NaiveBayes.cs | 8 +- .../Algorithms/Tabular/SARSA.cs | 3 +- docs/AdvancedGameEngineRoadMap.md | 218 ------------ docs/ExoplanetEngineRoadMap.md | 244 -------------- docs/Multiphysics-Roadmap.md | 140 -------- docs/Roadmap-v4.3.md | 16 +- docs/TerrainSpreadRoadMap.md | 319 ------------------ docs/completed/Multiphysics-Roadmap.md | 14 +- 10 files changed, 138 insertions(+), 941 deletions(-) create mode 100644 Numerics/NumericTest/NaiveBayesTests.cs delete mode 100644 docs/AdvancedGameEngineRoadMap.md delete mode 100644 docs/ExoplanetEngineRoadMap.md delete mode 100644 docs/Multiphysics-Roadmap.md delete mode 100644 docs/TerrainSpreadRoadMap.md 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/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/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/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 6e55216..feedaf3 100644 --- a/docs/Roadmap-v4.3.md +++ b/docs/Roadmap-v4.3.md @@ -240,11 +240,17 @@ Punkter som stod i v4.1-scopet och ännu inte är gjorda: 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 -- [ ] `NaiveBayes.NumClasses` sätts i `Fit` -- [ ] Ta bort tempfilsreferensen i csproj -- [ ] Ta bort dubblerade roadmaps ur `docs/` respektive `docs/completed/` -- [ ] Ta bort `SARSA._pendingNextAction` +### 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/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