From 51e529e07eaf135972c008a088b6cf8a9f755eca Mon Sep 17 00:00:00 2001 From: Candy Chiu Date: Tue, 30 Apr 2013 17:07:08 -0400 Subject: [PATCH 1/9] added BrentRootFinder --- src/Numerics/RootFinders/BrentRootFinder.cs | 114 +++++++++++++ src/Numerics/RootFinders/RootFinder.cs | 175 ++++++++++++++++++++ src/Numerics/RootFinders/RootFinderTest.cs | 24 +++ 3 files changed, 313 insertions(+) create mode 100644 src/Numerics/RootFinders/BrentRootFinder.cs create mode 100644 src/Numerics/RootFinders/RootFinder.cs create mode 100644 src/Numerics/RootFinders/RootFinderTest.cs diff --git a/src/Numerics/RootFinders/BrentRootFinder.cs b/src/Numerics/RootFinders/BrentRootFinder.cs new file mode 100644 index 00000000..06c926c6 --- /dev/null +++ b/src/Numerics/RootFinders/BrentRootFinder.cs @@ -0,0 +1,114 @@ +using System; + +namespace MathNet.Numerics.RootFinders +{ + public class BrentRootFinder : RootFinder + { + public BrentRootFinder() : base() + { + } + public BrentRootFinder(int numIters, double accuracy) : base(numIters, accuracy) + { + } + + protected override double Find() + { + /* The implementation of the algorithm was inspired by + Press, Teukolsky, Vetterling, and Flannery, + "Numerical Recipes in C", 2nd edition, Cambridge + University Press + */ + + double min1, min2; + double p, q, r, s, xAcc1, xMid = 0; + double d = 0.0, e = 0.0; + + // set up + double xmin = XMin; + double fxmin = Func(XMin); + double xmax = XMax; + double fxmax = Func(XMax); + + double root = xmax; + double froot = fxmax; + + // solve + int i = 0; + for (; i <= Iterations; i++) + { + if (Math.Sign(froot) == Math.Sign(fxmax)) + { + // Rename xMin_, root_, xMax_ and adjust bounds + xmax = xmin; + fxmax = fxmin; + e = d = root - xmin; + } + if (Math.Abs(fxmax) < Math.Abs(froot)) + { + xmin = root; + root = xmax; + xmax = xmin; + fxmin = froot; + froot = fxmax; + fxmax = fxmin; + } + // Convergence check + xAcc1 = 2.0 * DOUBLE_ACCURACY * Math.Abs(root) + 0.5 * Accuracy; + xMid = (xmax - root) / 2.0; + if (Math.Abs(xMid) <= xAcc1 || Close(froot, 0.0)) + { + return root; + } + if (Math.Abs(e) >= xAcc1 && + Math.Abs(fxmin) > Math.Abs(froot)) + { + + // Attempt inverse quadratic interpolation + s = froot / fxmin; + if (Close(xmin, xmax)) + { + p = 2.0 * xMid * s; + q = 1.0 - s; + } + else + { + q = fxmin / fxmax; + r = froot / fxmax; + p = s * (2.0 * xMid * q * (q - r) - (root - xmin) * (r - 1.0)); + q = (q - 1.0) * (r - 1.0) * (s - 1.0); + } + if (p > 0.0) q = -q; // Check whether in bounds + p = Math.Abs(p); + min1 = 3.0 * xMid * q - Math.Abs(xAcc1 * q); + min2 = Math.Abs(e * q); + if (2.0 * p < Math.Min(min1, min2)) + { + e = d; // Accept interpolation + d = p / q; + } + else + { + d = xMid; // Interpolation failed, use bisection + e = d; + } + } + else + { + // Bounds decreasing too slowly, use bisection + d = xMid; + e = d; + } + xmin = root; + fxmin = froot; + if (Math.Abs(d) > xAcc1) + root += d; + else + root += Sign(xAcc1, xMid); + froot = Func(root); + } + + // The algorithm has exceeded the number of iterations allowed + throw new RootFinderException(ACCURACY_NOT_REACHED, i, new Range(XMin, XMax), Math.Abs(xMid)); + } + } +} diff --git a/src/Numerics/RootFinders/RootFinder.cs b/src/Numerics/RootFinders/RootFinder.cs new file mode 100644 index 00000000..4ee24498 --- /dev/null +++ b/src/Numerics/RootFinders/RootFinder.cs @@ -0,0 +1,175 @@ +using System; + +namespace MathNet.Numerics.RootFinders +{ + public struct Range + { + double Min, Max; + + public Range(double min, double max) + { + Min = min; Max = max; + } + } + public class RootFinderException : Exception + { + private int m_Iteration; + private Range m_Range; + private double m_Accuracy; + + public RootFinderException(string message, int iteration, Range range, double accuracy) + : base(message) + { + m_Iteration = iteration; + m_Range = range; + m_Accuracy = accuracy; + } + + public int Iteration + { + get { return m_Iteration; } + set { m_Iteration = value; } + } + + public Range Range + { + get { return m_Range; } + set { m_Range = value; } + } + + public double Accuracy { set; get; } + } + public abstract class RootFinder + { + + protected const string INVALID_RANGE="Invalid range while finding root"; + protected const string ACCURACY_NOT_REACHED = "The accuracy couldn't be reached with the specified number of iterations"; + protected const string ROOT_NOT_FOUND = "The algorithm ended without root in the range"; + protected const string ROOT_NOT_BRACKETED = "The algorithm could not start because the root seemed not to be bracketed"; + protected const string INVALID_ALGORITHM = "This algorithm is not able to solve this equation"; + protected const double DOUBLE_ACCURACY = 9.99200722162641E-16; + private const int DEFAULT_MAX_ITERATIONS = 30; + private const double DEFAULT_ACCURACY = 1e-8; + + //protected int _maxNumIters; + //protected double _xmin = double.MinValue; + //protected double _xmax = double.MaxValue; + //protected double _accuracy; + //protected Func _func; + //protected Func m_Of; + int _maxNumIters; + double _xmin = double.MinValue; + double _xmax = double.MaxValue; + double _accuracy; + Func _func; + private double bracketingFactor = 1.6; + + + + /// Constructor. + /// A continuous function. + public RootFinder() : this(DEFAULT_MAX_ITERATIONS, DEFAULT_ACCURACY) + { + } + + public RootFinder(int numIters, double accuracy) + { + _maxNumIters = numIters; + _accuracy = accuracy; + } + + #region Properties + protected double XMin { get { return _xmin; } } + protected double XMax { get { return _xmax; } } + + public Func Func + { + get { return _func; } + set { _func = value; } + } + public double BracketingFactor + { + get { return bracketingFactor; } + set + { + if (value <= 0.0) throw new ArgumentOutOfRangeException(); + bracketingFactor = value; + } + } + public int Iterations + { + set + { + if (value <= 0) throw new ArgumentOutOfRangeException(); + _maxNumIters = value; + } + protected get { return _maxNumIters; } + } + public double Accuracy + { + get { return _accuracy; } + set { _accuracy = value; } + } + #endregion Properties + + /// Detect a range containing at least one root. + /// Lower value of the range. + /// Upper value of the range + /// The growing factor of research. Usually 1.6. + /// True if the bracketing operation succeeded, else otherwise. + /// This iterative methods stops when two values with opposite signs are found. + public bool SearchBracketsOutward(ref double xmin, ref double xmax, double factor) + { + if (xmin >= xmax) + { + throw new RootFinderException(INVALID_RANGE, 0, new Range(xmin, xmax), 0.0); + } + + double fmin = _func(xmin); + double fmax = _func(xmax); + + int i = 0; + while (i++ < _maxNumIters) + { + if (Math.Sign(fmin) != Math.Sign(fmax)) return true; + if (Math.Abs(fmin) < Math.Abs(fmax)) + { + xmin += factor * (xmin - xmax); + fmin = _func(xmin); + } + else + { + xmax += factor * (xmax - xmin); + fmax = _func(xmax); + } + } + + throw new RootFinderException(ROOT_NOT_FOUND, i, new Range(fmin, fmax), 0.0); + } + + /// Prototype algorithm for solving the equation f(x)=0. + /// The low value of the range where the root is supposed to be. + /// The high value of the range where the root is supposed to be. + /// Returns the root with the specified accuracy. + public virtual double Solve(double x1, double x2) + { + _xmin = x1; + _xmax = x2; + return Find(); + } + + protected abstract double Find(); + + /// Helper method useful for preventing rounding errors. + /// a*sign(b) + protected static double Sign(double a, double b) + { + return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); + } + + protected static bool Close(double d1, double d2) + { + return Math.Abs(d1 - d2) <= double.Epsilon; + } + } +} diff --git a/src/Numerics/RootFinders/RootFinderTest.cs b/src/Numerics/RootFinders/RootFinderTest.cs new file mode 100644 index 00000000..02f731de --- /dev/null +++ b/src/Numerics/RootFinders/RootFinderTest.cs @@ -0,0 +1,24 @@ +namespace MathNet.Numerics.RootFinders +{ + using System; + using System.Collections.Generic; + using MbUnit.Framework; + using Gallio.Framework; + using SuntrustPortfolio.Numerics; + + [TestFixture] + public class RootFinderTest + { + BrentRootFinder _solver = new BrentRootFinder(100, 1e-14); + + [Test] + public void MultipleRoots() + { + Func f = (x) => { return x * x - 4; }; + _solver.Func = f; + double root = _solver.Solve(-5, 5); + + Assert.AreEqual(0, f(root)); + } + } +} From be5d8b9ebc914a5d852ae152f0c0e1703e5c08d0 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Wed, 1 May 2013 17:29:30 +0200 Subject: [PATCH 2/9] RootFinding: move files, namespace, tests, projects --- src/Numerics/Numerics.csproj | 2 ++ .../BrentRootFinder.cs | 2 +- .../{RootFinders => RootFinding}/RootFinder.cs | 2 +- .../RootFindingTests/BrentRootFinderTest.cs} | 16 +++++++--------- src/UnitTests/UnitTests.csproj | 1 + 5 files changed, 12 insertions(+), 11 deletions(-) rename src/Numerics/{RootFinders => RootFinding}/BrentRootFinder.cs (98%) rename src/Numerics/{RootFinders => RootFinding}/RootFinder.cs (99%) rename src/{Numerics/RootFinders/RootFinderTest.cs => UnitTests/RootFindingTests/BrentRootFinderTest.cs} (50%) diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 9249bd1a..b44a866c 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -109,6 +109,8 @@ + + diff --git a/src/Numerics/RootFinders/BrentRootFinder.cs b/src/Numerics/RootFinding/BrentRootFinder.cs similarity index 98% rename from src/Numerics/RootFinders/BrentRootFinder.cs rename to src/Numerics/RootFinding/BrentRootFinder.cs index 06c926c6..82ec8b2d 100644 --- a/src/Numerics/RootFinders/BrentRootFinder.cs +++ b/src/Numerics/RootFinding/BrentRootFinder.cs @@ -1,6 +1,6 @@ using System; -namespace MathNet.Numerics.RootFinders +namespace MathNet.Numerics.RootFinding { public class BrentRootFinder : RootFinder { diff --git a/src/Numerics/RootFinders/RootFinder.cs b/src/Numerics/RootFinding/RootFinder.cs similarity index 99% rename from src/Numerics/RootFinders/RootFinder.cs rename to src/Numerics/RootFinding/RootFinder.cs index 4ee24498..94115c6f 100644 --- a/src/Numerics/RootFinders/RootFinder.cs +++ b/src/Numerics/RootFinding/RootFinder.cs @@ -1,6 +1,6 @@ using System; -namespace MathNet.Numerics.RootFinders +namespace MathNet.Numerics.RootFinding { public struct Range { diff --git a/src/Numerics/RootFinders/RootFinderTest.cs b/src/UnitTests/RootFindingTests/BrentRootFinderTest.cs similarity index 50% rename from src/Numerics/RootFinders/RootFinderTest.cs rename to src/UnitTests/RootFindingTests/BrentRootFinderTest.cs index 02f731de..9c00ee76 100644 --- a/src/Numerics/RootFinders/RootFinderTest.cs +++ b/src/UnitTests/RootFindingTests/BrentRootFinderTest.cs @@ -1,15 +1,13 @@ -namespace MathNet.Numerics.RootFinders -{ - using System; - using System.Collections.Generic; - using MbUnit.Framework; - using Gallio.Framework; - using SuntrustPortfolio.Numerics; +using MathNet.Numerics.RootFinding; +using NUnit.Framework; +using System; +namespace MathNet.Numerics.UnitTests.RootFindingTests +{ [TestFixture] - public class RootFinderTest + public class BrentRootFinderTest { - BrentRootFinder _solver = new BrentRootFinder(100, 1e-14); + readonly BrentRootFinder _solver = new BrentRootFinder(100, 1e-14); [Test] public void MultipleRoots() diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 1307db01..552e2bef 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -770,6 +770,7 @@ + From 8742c9257835286ca7d7e87ba46978262f9909b3 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Wed, 1 May 2013 17:39:44 +0200 Subject: [PATCH 3/9] RootFinding: migrate message strings to resources --- src/Numerics/Properties/Resources.Designer.cs | 18 ++++++++++++++++++ src/Numerics/Properties/Resources.resx | 14 ++++++++++---- src/Numerics/RootFinding/BrentRootFinder.cs | 3 ++- src/Numerics/RootFinding/RootFinder.cs | 13 ++++--------- 4 files changed, 34 insertions(+), 14 deletions(-) diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs index 13e9f0d2..0a336eb0 100644 --- a/src/Numerics/Properties/Resources.Designer.cs +++ b/src/Numerics/Properties/Resources.Designer.cs @@ -60,6 +60,15 @@ namespace MathNet.Numerics.Properties { } } + /// + /// Looks up a localized string similar to The accuracy couldn't be reached with the specified number of iterations.. + /// + public static string AccuracyNotReached { + get { + return ResourceManager.GetString("AccuracyNotReached", resourceCulture); + } + } + /// /// Looks up a localized string similar to The array arguments must have the same length.. /// @@ -744,6 +753,15 @@ namespace MathNet.Numerics.Properties { } } + /// + /// Looks up a localized string similar to The algorithm ended without root in the range.. + /// + public static string RootNotFound { + get { + return ResourceManager.GetString("RootNotFound", resourceCulture); + } + } + /// /// Looks up a localized string similar to The number of rows must greater than or equal to the number of columns.. /// diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx index 2d7c5c10..e0645a91 100644 --- a/src/Numerics/Properties/Resources.resx +++ b/src/Numerics/Properties/Resources.resx @@ -369,10 +369,16 @@ We only support sparse matrix with less than int.MaxValue elements. - - Sample points should be sorted in strictly ascending order + + Sample points should be sorted in strictly ascending order - - All sample points should be unique. + + All sample points should be unique. + + + The accuracy couldn't be reached with the specified number of iterations. + + + The algorithm ended without root in the range. \ No newline at end of file diff --git a/src/Numerics/RootFinding/BrentRootFinder.cs b/src/Numerics/RootFinding/BrentRootFinder.cs index 82ec8b2d..d5f4bb90 100644 --- a/src/Numerics/RootFinding/BrentRootFinder.cs +++ b/src/Numerics/RootFinding/BrentRootFinder.cs @@ -1,4 +1,5 @@ using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.RootFinding { @@ -108,7 +109,7 @@ namespace MathNet.Numerics.RootFinding } // The algorithm has exceeded the number of iterations allowed - throw new RootFinderException(ACCURACY_NOT_REACHED, i, new Range(XMin, XMax), Math.Abs(xMid)); + throw new RootFinderException(Resources.AccuracyNotReached, i, new Range(XMin, XMax), Math.Abs(xMid)); } } } diff --git a/src/Numerics/RootFinding/RootFinder.cs b/src/Numerics/RootFinding/RootFinder.cs index 94115c6f..4d56a650 100644 --- a/src/Numerics/RootFinding/RootFinder.cs +++ b/src/Numerics/RootFinding/RootFinder.cs @@ -1,4 +1,5 @@ using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.RootFinding { @@ -41,13 +42,7 @@ namespace MathNet.Numerics.RootFinding } public abstract class RootFinder { - - protected const string INVALID_RANGE="Invalid range while finding root"; - protected const string ACCURACY_NOT_REACHED = "The accuracy couldn't be reached with the specified number of iterations"; - protected const string ROOT_NOT_FOUND = "The algorithm ended without root in the range"; - protected const string ROOT_NOT_BRACKETED = "The algorithm could not start because the root seemed not to be bracketed"; - protected const string INVALID_ALGORITHM = "This algorithm is not able to solve this equation"; - protected const double DOUBLE_ACCURACY = 9.99200722162641E-16; + protected const double DOUBLE_ACCURACY = 9.99200722162641E-16; private const int DEFAULT_MAX_ITERATIONS = 30; private const double DEFAULT_ACCURACY = 1e-8; @@ -122,7 +117,7 @@ namespace MathNet.Numerics.RootFinding { if (xmin >= xmax) { - throw new RootFinderException(INVALID_RANGE, 0, new Range(xmin, xmax), 0.0); + throw new RootFinderException(string.Format(Resources.ArgumentOutOfRangeGreater,"xmax","xmin"), 0, new Range(xmin, xmax), 0.0); } double fmin = _func(xmin); @@ -144,7 +139,7 @@ namespace MathNet.Numerics.RootFinding } } - throw new RootFinderException(ROOT_NOT_FOUND, i, new Range(fmin, fmax), 0.0); + throw new RootFinderException(Resources.RootNotFound, i, new Range(fmin, fmax), 0.0); } /// Prototype algorithm for solving the equation f(x)=0. From 2d4ac98a2971c0628d476bfea4f4a9cbb311aa9f Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Wed, 1 May 2013 17:44:31 +0200 Subject: [PATCH 4/9] RootFinding: simplify exception --- src/Numerics/Numerics.csproj | 1 + src/Numerics/RootFinding/BrentRootFinder.cs | 2 +- src/Numerics/RootFinding/RootFinder.cs | 45 ++----------------- .../RootFinding/RootFindingException.cs | 21 +++++++++ 4 files changed, 26 insertions(+), 43 deletions(-) create mode 100644 src/Numerics/RootFinding/RootFindingException.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index b44a866c..23c22630 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -111,6 +111,7 @@ + diff --git a/src/Numerics/RootFinding/BrentRootFinder.cs b/src/Numerics/RootFinding/BrentRootFinder.cs index d5f4bb90..364cf333 100644 --- a/src/Numerics/RootFinding/BrentRootFinder.cs +++ b/src/Numerics/RootFinding/BrentRootFinder.cs @@ -109,7 +109,7 @@ namespace MathNet.Numerics.RootFinding } // The algorithm has exceeded the number of iterations allowed - throw new RootFinderException(Resources.AccuracyNotReached, i, new Range(XMin, XMax), Math.Abs(xMid)); + throw new RootFindingException(Resources.AccuracyNotReached, i, XMin, XMax, Math.Abs(xMid)); } } } diff --git a/src/Numerics/RootFinding/RootFinder.cs b/src/Numerics/RootFinding/RootFinder.cs index 4d56a650..8fe6540d 100644 --- a/src/Numerics/RootFinding/RootFinder.cs +++ b/src/Numerics/RootFinding/RootFinder.cs @@ -3,44 +3,7 @@ using MathNet.Numerics.Properties; namespace MathNet.Numerics.RootFinding { - public struct Range - { - double Min, Max; - - public Range(double min, double max) - { - Min = min; Max = max; - } - } - public class RootFinderException : Exception - { - private int m_Iteration; - private Range m_Range; - private double m_Accuracy; - - public RootFinderException(string message, int iteration, Range range, double accuracy) - : base(message) - { - m_Iteration = iteration; - m_Range = range; - m_Accuracy = accuracy; - } - - public int Iteration - { - get { return m_Iteration; } - set { m_Iteration = value; } - } - - public Range Range - { - get { return m_Range; } - set { m_Range = value; } - } - - public double Accuracy { set; get; } - } - public abstract class RootFinder + public abstract class RootFinder { protected const double DOUBLE_ACCURACY = 9.99200722162641E-16; private const int DEFAULT_MAX_ITERATIONS = 30; @@ -59,8 +22,6 @@ namespace MathNet.Numerics.RootFinding Func _func; private double bracketingFactor = 1.6; - - /// Constructor. /// A continuous function. public RootFinder() : this(DEFAULT_MAX_ITERATIONS, DEFAULT_ACCURACY) @@ -117,7 +78,7 @@ namespace MathNet.Numerics.RootFinding { if (xmin >= xmax) { - throw new RootFinderException(string.Format(Resources.ArgumentOutOfRangeGreater,"xmax","xmin"), 0, new Range(xmin, xmax), 0.0); + throw new RootFindingException(string.Format(Resources.ArgumentOutOfRangeGreater,"xmax","xmin"), 0, xmin, xmax, 0.0); } double fmin = _func(xmin); @@ -139,7 +100,7 @@ namespace MathNet.Numerics.RootFinding } } - throw new RootFinderException(Resources.RootNotFound, i, new Range(fmin, fmax), 0.0); + throw new RootFindingException(Resources.RootNotFound, i, fmin, fmax, 0.0); } /// Prototype algorithm for solving the equation f(x)=0. diff --git a/src/Numerics/RootFinding/RootFindingException.cs b/src/Numerics/RootFinding/RootFindingException.cs new file mode 100644 index 00000000..3aaac181 --- /dev/null +++ b/src/Numerics/RootFinding/RootFindingException.cs @@ -0,0 +1,21 @@ +using System; + +namespace MathNet.Numerics.RootFinding +{ + public class RootFindingException : Exception + { + public RootFindingException(string message, int iteration, double rangeMin, double rangeMax, double accuracy) + : base(message) + { + Iteration = iteration; + RangeMin = rangeMin; + RangeMax = rangeMax; + Accuracy = accuracy; + } + + public int Iteration { get; set; } + public double RangeMin { get; set; } + public double RangeMax { get; set; } + public double Accuracy { set; get; } + } +} \ No newline at end of file From a177228238eb5113fa04dcd7d70e513ece402808 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Wed, 1 May 2013 17:47:14 +0200 Subject: [PATCH 5/9] RootFinding: simplify unit test --- .../RootFindingTests/BrentRootFinderTest.cs | 22 ------------------- .../RootFindingTests/BrentRootFindingTest.cs | 17 ++++++++++++++ src/UnitTests/UnitTests.csproj | 2 +- 3 files changed, 18 insertions(+), 23 deletions(-) delete mode 100644 src/UnitTests/RootFindingTests/BrentRootFinderTest.cs create mode 100644 src/UnitTests/RootFindingTests/BrentRootFindingTest.cs diff --git a/src/UnitTests/RootFindingTests/BrentRootFinderTest.cs b/src/UnitTests/RootFindingTests/BrentRootFinderTest.cs deleted file mode 100644 index 9c00ee76..00000000 --- a/src/UnitTests/RootFindingTests/BrentRootFinderTest.cs +++ /dev/null @@ -1,22 +0,0 @@ -using MathNet.Numerics.RootFinding; -using NUnit.Framework; -using System; - -namespace MathNet.Numerics.UnitTests.RootFindingTests -{ - [TestFixture] - public class BrentRootFinderTest - { - readonly BrentRootFinder _solver = new BrentRootFinder(100, 1e-14); - - [Test] - public void MultipleRoots() - { - Func f = (x) => { return x * x - 4; }; - _solver.Func = f; - double root = _solver.Solve(-5, 5); - - Assert.AreEqual(0, f(root)); - } - } -} diff --git a/src/UnitTests/RootFindingTests/BrentRootFindingTest.cs b/src/UnitTests/RootFindingTests/BrentRootFindingTest.cs new file mode 100644 index 00000000..4a264a91 --- /dev/null +++ b/src/UnitTests/RootFindingTests/BrentRootFindingTest.cs @@ -0,0 +1,17 @@ +using MathNet.Numerics.RootFinding; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.RootFindingTests +{ + [TestFixture] + public class BrentRootFindingTest + { + [Test] + public void MultipleRoots() + { + var solver = new BrentRootFinder(100, 1e-14) {Func = x => x*x - 4}; + double root = solver.Solve(-5, 5); + Assert.AreEqual(0, solver.Func(root)); + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 552e2bef..916dc10f 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -770,7 +770,7 @@ - + From 65cfd6c46eb8a2bb922cae76fb68166036c23770 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Wed, 1 May 2013 17:50:23 +0200 Subject: [PATCH 6/9] RootFinding: cosmetics --- src/Numerics/RootFinding/BrentRootFinder.cs | 2 +- src/Numerics/RootFinding/RootFinder.cs | 38 +++++++-------------- 2 files changed, 14 insertions(+), 26 deletions(-) diff --git a/src/Numerics/RootFinding/BrentRootFinder.cs b/src/Numerics/RootFinding/BrentRootFinder.cs index 364cf333..2241a61f 100644 --- a/src/Numerics/RootFinding/BrentRootFinder.cs +++ b/src/Numerics/RootFinding/BrentRootFinder.cs @@ -54,7 +54,7 @@ namespace MathNet.Numerics.RootFinding fxmax = fxmin; } // Convergence check - xAcc1 = 2.0 * DOUBLE_ACCURACY * Math.Abs(root) + 0.5 * Accuracy; + xAcc1 = 2.0 * DoubleAccuracy * Math.Abs(root) + 0.5 * Accuracy; xMid = (xmax - root) / 2.0; if (Math.Abs(xMid) <= xAcc1 || Close(froot, 0.0)) { diff --git a/src/Numerics/RootFinding/RootFinder.cs b/src/Numerics/RootFinding/RootFinder.cs index 8fe6540d..024bcac4 100644 --- a/src/Numerics/RootFinding/RootFinder.cs +++ b/src/Numerics/RootFinding/RootFinder.cs @@ -5,36 +5,26 @@ namespace MathNet.Numerics.RootFinding { public abstract class RootFinder { - protected const double DOUBLE_ACCURACY = 9.99200722162641E-16; - private const int DEFAULT_MAX_ITERATIONS = 30; - private const double DEFAULT_ACCURACY = 1e-8; - - //protected int _maxNumIters; - //protected double _xmin = double.MinValue; - //protected double _xmax = double.MaxValue; - //protected double _accuracy; - //protected Func _func; - //protected Func m_Of; + protected const double DoubleAccuracy = 9.99200722162641E-16; + private const int DefaultMaxIterations = 30; + private const double DefaultAccuracy = 1e-8; + int _maxNumIters; double _xmin = double.MinValue; double _xmax = double.MaxValue; - double _accuracy; Func _func; - private double bracketingFactor = 1.6; + private double _bracketingFactor = 1.6; - /// Constructor. - /// A continuous function. - public RootFinder() : this(DEFAULT_MAX_ITERATIONS, DEFAULT_ACCURACY) + public RootFinder() : this(DefaultMaxIterations, DefaultAccuracy) { } public RootFinder(int numIters, double accuracy) { _maxNumIters = numIters; - _accuracy = accuracy; + Accuracy = accuracy; } - #region Properties protected double XMin { get { return _xmin; } } protected double XMax { get { return _xmax; } } @@ -43,15 +33,19 @@ namespace MathNet.Numerics.RootFinding get { return _func; } set { _func = value; } } + + public double Accuracy { get; set; } + public double BracketingFactor { - get { return bracketingFactor; } + get { return _bracketingFactor; } set { if (value <= 0.0) throw new ArgumentOutOfRangeException(); - bracketingFactor = value; + _bracketingFactor = value; } } + public int Iterations { set @@ -61,12 +55,6 @@ namespace MathNet.Numerics.RootFinding } protected get { return _maxNumIters; } } - public double Accuracy - { - get { return _accuracy; } - set { _accuracy = value; } - } - #endregion Properties /// Detect a range containing at least one root. /// Lower value of the range. From 2f9913921cd27be3a9472cd528f416f098e2f096 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Wed, 1 May 2013 18:30:00 +0200 Subject: [PATCH 7/9] RootFinding: extract bracketing to separate class --- src/Numerics/Numerics.csproj | 1 + src/Numerics/RootFinding/Bracketing.cs | 48 +++++++++++++++++ src/Numerics/RootFinding/BrentRootFinder.cs | 12 +++++ src/Numerics/RootFinding/RootFinder.cs | 59 --------------------- 4 files changed, 61 insertions(+), 59 deletions(-) create mode 100644 src/Numerics/RootFinding/Bracketing.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 23c22630..f3f830c6 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -109,6 +109,7 @@ + diff --git a/src/Numerics/RootFinding/Bracketing.cs b/src/Numerics/RootFinding/Bracketing.cs new file mode 100644 index 00000000..0cc1f603 --- /dev/null +++ b/src/Numerics/RootFinding/Bracketing.cs @@ -0,0 +1,48 @@ +using System; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.RootFinding +{ + public static class Bracketing + { + /// Detect a range containing at least one root. + /// The function to detect roots from. + /// Lower value of the range. + /// Upper value of the range + /// The growing factor of research. Usually 1.6. + /// Maximum number of iterations. Usually 50. + /// True if the bracketing operation succeeded, false otherwise. + /// This iterative methods stops when two values with opposite signs are found. + public static bool SearchOutward(Func f, ref double xmin, ref double xmax, double factor = 1.6, int maxIterations = 50) + { + if (xmin >= xmax) + { + throw new ArgumentOutOfRangeException("xmax", string.Format(Resources.ArgumentOutOfRangeGreater, "xmax", "xmin")); + } + + double fmin = f(xmin); + double fmax = f(xmax); + + for(int i=0;iHelper method useful for preventing rounding errors. + /// a*sign(b) + static double Sign(double a, double b) + { + return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); + } + + static bool Close(double d1, double d2) + { + return Math.Abs(d1 - d2) <= double.Epsilon; + } } } diff --git a/src/Numerics/RootFinding/RootFinder.cs b/src/Numerics/RootFinding/RootFinder.cs index 024bcac4..71c3e8f0 100644 --- a/src/Numerics/RootFinding/RootFinder.cs +++ b/src/Numerics/RootFinding/RootFinder.cs @@ -1,5 +1,4 @@ using System; -using MathNet.Numerics.Properties; namespace MathNet.Numerics.RootFinding { @@ -13,7 +12,6 @@ namespace MathNet.Numerics.RootFinding double _xmin = double.MinValue; double _xmax = double.MaxValue; Func _func; - private double _bracketingFactor = 1.6; public RootFinder() : this(DefaultMaxIterations, DefaultAccuracy) { @@ -36,16 +34,6 @@ namespace MathNet.Numerics.RootFinding public double Accuracy { get; set; } - public double BracketingFactor - { - get { return _bracketingFactor; } - set - { - if (value <= 0.0) throw new ArgumentOutOfRangeException(); - _bracketingFactor = value; - } - } - public int Iterations { set @@ -56,41 +44,6 @@ namespace MathNet.Numerics.RootFinding protected get { return _maxNumIters; } } - /// Detect a range containing at least one root. - /// Lower value of the range. - /// Upper value of the range - /// The growing factor of research. Usually 1.6. - /// True if the bracketing operation succeeded, else otherwise. - /// This iterative methods stops when two values with opposite signs are found. - public bool SearchBracketsOutward(ref double xmin, ref double xmax, double factor) - { - if (xmin >= xmax) - { - throw new RootFindingException(string.Format(Resources.ArgumentOutOfRangeGreater,"xmax","xmin"), 0, xmin, xmax, 0.0); - } - - double fmin = _func(xmin); - double fmax = _func(xmax); - - int i = 0; - while (i++ < _maxNumIters) - { - if (Math.Sign(fmin) != Math.Sign(fmax)) return true; - if (Math.Abs(fmin) < Math.Abs(fmax)) - { - xmin += factor * (xmin - xmax); - fmin = _func(xmin); - } - else - { - xmax += factor * (xmax - xmin); - fmax = _func(xmax); - } - } - - throw new RootFindingException(Resources.RootNotFound, i, fmin, fmax, 0.0); - } - /// Prototype algorithm for solving the equation f(x)=0. /// The low value of the range where the root is supposed to be. /// The high value of the range where the root is supposed to be. @@ -103,17 +56,5 @@ namespace MathNet.Numerics.RootFinding } protected abstract double Find(); - - /// Helper method useful for preventing rounding errors. - /// a*sign(b) - protected static double Sign(double a, double b) - { - return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); - } - - protected static bool Close(double d1, double d2) - { - return Math.Abs(d1 - d2) <= double.Epsilon; - } } } From 04a5f30d99014f104bcea27ae2e2bd6cdf4dda2e Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Wed, 1 May 2013 18:49:24 +0200 Subject: [PATCH 8/9] RootFinding: refactoring, simplify --- src/Numerics/Numerics.csproj | 3 +- .../{BrentRootFinder.cs => FindRoots.cs} | 45 ++++++-------- src/Numerics/RootFinding/RootFinder.cs | 60 ------------------- .../{BrentRootFindingTest.cs => BrentTest.cs} | 7 +-- src/UnitTests/UnitTests.csproj | 2 +- 5 files changed, 24 insertions(+), 93 deletions(-) rename src/Numerics/RootFinding/{BrentRootFinder.cs => FindRoots.cs} (70%) delete mode 100644 src/Numerics/RootFinding/RootFinder.cs rename src/UnitTests/RootFindingTests/{BrentRootFindingTest.cs => BrentTest.cs} (50%) diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index f3f830c6..2eee6c0b 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -110,8 +110,7 @@ - - + diff --git a/src/Numerics/RootFinding/BrentRootFinder.cs b/src/Numerics/RootFinding/FindRoots.cs similarity index 70% rename from src/Numerics/RootFinding/BrentRootFinder.cs rename to src/Numerics/RootFinding/FindRoots.cs index 6f5b9096..adf93d84 100644 --- a/src/Numerics/RootFinding/BrentRootFinder.cs +++ b/src/Numerics/RootFinding/FindRoots.cs @@ -3,39 +3,32 @@ using MathNet.Numerics.Properties; namespace MathNet.Numerics.RootFinding { - public class BrentRootFinder : RootFinder + public static class FindRoots { - public BrentRootFinder() : base() + /// Find a solution of the equation f(x)=0. + /// The function to find roots from. + /// The low value of the range where the root is supposed to be. + /// The high value of the range where the root is supposed to be. + /// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. + /// Maximum number of iterations. Usually 100. + /// Returns the root with the specified accuracy. + /// + /// Algorithm by by Brent, Van Wijngaarden, Dekker et al. + /// Implementation inspired by Press, Teukolsky, Vetterling, and Flannery, "Numerical Recipes in C", 2nd edition, Cambridge University Press + /// + public static double BrentMethod(Func f, double xmin, double xmax, double accuracy = 1e-8, int maxIterations = 100) { - } - public BrentRootFinder(int numIters, double accuracy) : base(numIters, accuracy) - { - } - - protected override double Find() - { - /* The implementation of the algorithm was inspired by - Press, Teukolsky, Vetterling, and Flannery, - "Numerical Recipes in C", 2nd edition, Cambridge - University Press - */ - double min1, min2; double p, q, r, s, xAcc1, xMid = 0; double d = 0.0, e = 0.0; // set up - double xmin = XMin; - double fxmin = Func(XMin); - double xmax = XMax; - double fxmax = Func(XMax); - + double fxmin = f(xmin); + double fxmax = f(xmax); double root = xmax; double froot = fxmax; - // solve - int i = 0; - for (; i <= Iterations; i++) + for (int i = 0; i <= maxIterations; i++) { if (Math.Sign(froot) == Math.Sign(fxmax)) { @@ -54,7 +47,7 @@ namespace MathNet.Numerics.RootFinding fxmax = fxmin; } // Convergence check - xAcc1 = 2.0 * DoubleAccuracy * Math.Abs(root) + 0.5 * Accuracy; + xAcc1 = 2.0 * Precision.DoubleMachinePrecision * Math.Abs(root) + 0.5 * accuracy; xMid = (xmax - root) / 2.0; if (Math.Abs(xMid) <= xAcc1 || Close(froot, 0.0)) { @@ -105,11 +98,11 @@ namespace MathNet.Numerics.RootFinding root += d; else root += Sign(xAcc1, xMid); - froot = Func(root); + froot = f(root); } // The algorithm has exceeded the number of iterations allowed - throw new RootFindingException(Resources.AccuracyNotReached, i, XMin, XMax, Math.Abs(xMid)); + throw new RootFindingException(Resources.AccuracyNotReached, maxIterations, xmin, xmax, Math.Abs(xMid)); } /// Helper method useful for preventing rounding errors. diff --git a/src/Numerics/RootFinding/RootFinder.cs b/src/Numerics/RootFinding/RootFinder.cs deleted file mode 100644 index 71c3e8f0..00000000 --- a/src/Numerics/RootFinding/RootFinder.cs +++ /dev/null @@ -1,60 +0,0 @@ -using System; - -namespace MathNet.Numerics.RootFinding -{ - public abstract class RootFinder - { - protected const double DoubleAccuracy = 9.99200722162641E-16; - private const int DefaultMaxIterations = 30; - private const double DefaultAccuracy = 1e-8; - - int _maxNumIters; - double _xmin = double.MinValue; - double _xmax = double.MaxValue; - Func _func; - - public RootFinder() : this(DefaultMaxIterations, DefaultAccuracy) - { - } - - public RootFinder(int numIters, double accuracy) - { - _maxNumIters = numIters; - Accuracy = accuracy; - } - - protected double XMin { get { return _xmin; } } - protected double XMax { get { return _xmax; } } - - public Func Func - { - get { return _func; } - set { _func = value; } - } - - public double Accuracy { get; set; } - - public int Iterations - { - set - { - if (value <= 0) throw new ArgumentOutOfRangeException(); - _maxNumIters = value; - } - protected get { return _maxNumIters; } - } - - /// Prototype algorithm for solving the equation f(x)=0. - /// The low value of the range where the root is supposed to be. - /// The high value of the range where the root is supposed to be. - /// Returns the root with the specified accuracy. - public virtual double Solve(double x1, double x2) - { - _xmin = x1; - _xmax = x2; - return Find(); - } - - protected abstract double Find(); - } -} diff --git a/src/UnitTests/RootFindingTests/BrentRootFindingTest.cs b/src/UnitTests/RootFindingTests/BrentTest.cs similarity index 50% rename from src/UnitTests/RootFindingTests/BrentRootFindingTest.cs rename to src/UnitTests/RootFindingTests/BrentTest.cs index 4a264a91..2956caad 100644 --- a/src/UnitTests/RootFindingTests/BrentRootFindingTest.cs +++ b/src/UnitTests/RootFindingTests/BrentTest.cs @@ -4,14 +4,13 @@ using NUnit.Framework; namespace MathNet.Numerics.UnitTests.RootFindingTests { [TestFixture] - public class BrentRootFindingTest + public class BrentTest { [Test] public void MultipleRoots() { - var solver = new BrentRootFinder(100, 1e-14) {Func = x => x*x - 4}; - double root = solver.Solve(-5, 5); - Assert.AreEqual(0, solver.Func(root)); + double root = FindRoots.BrentMethod(x => x*x - 4, -5, 5, 1e-14, 100); + Assert.AreEqual(0, root*root - 4); } } } diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 916dc10f..7c771f55 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -770,7 +770,7 @@ - + From ddac0a0b615b62caf276ee029759424aa8af8646 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Wed, 1 May 2013 19:26:24 +0200 Subject: [PATCH 9/9] RootFinding: algorithm cosmetics --- src/Numerics/RootFinding/FindRoots.cs | 55 +++++++++++++++------------ 1 file changed, 30 insertions(+), 25 deletions(-) diff --git a/src/Numerics/RootFinding/FindRoots.cs b/src/Numerics/RootFinding/FindRoots.cs index adf93d84..5604cabf 100644 --- a/src/Numerics/RootFinding/FindRoots.cs +++ b/src/Numerics/RootFinding/FindRoots.cs @@ -18,11 +18,9 @@ namespace MathNet.Numerics.RootFinding /// public static double BrentMethod(Func f, double xmin, double xmax, double accuracy = 1e-8, int maxIterations = 100) { - double min1, min2; - double p, q, r, s, xAcc1, xMid = 0; - double d = 0.0, e = 0.0; + double xMid = 0; + double d = 0.0, e = 0.0; - // set up double fxmin = f(xmin); double fxmax = f(xmax); double root = xmax; @@ -30,13 +28,14 @@ namespace MathNet.Numerics.RootFinding for (int i = 0; i <= maxIterations; i++) { + // adjust bounds if (Math.Sign(froot) == Math.Sign(fxmax)) { - // Rename xMin_, root_, xMax_ and adjust bounds xmax = xmin; fxmax = fxmin; e = d = root - xmin; } + if (Math.Abs(fxmax) < Math.Abs(froot)) { xmin = root; @@ -46,20 +45,22 @@ namespace MathNet.Numerics.RootFinding froot = fxmax; fxmax = fxmin; } - // Convergence check - xAcc1 = 2.0 * Precision.DoubleMachinePrecision * Math.Abs(root) + 0.5 * accuracy; + + // convergence check + double xAcc1 = 2.0 * Precision.DoubleMachinePrecision * Math.Abs(root) + 0.5 * accuracy; xMid = (xmax - root) / 2.0; - if (Math.Abs(xMid) <= xAcc1 || Close(froot, 0.0)) + if (Math.Abs(xMid) <= xAcc1 || froot.AlmostEqualWithAbsoluteError(0, froot, accuracy)) { return root; } - if (Math.Abs(e) >= xAcc1 && - Math.Abs(fxmin) > Math.Abs(froot)) - { + if (Math.Abs(e) >= xAcc1 && Math.Abs(fxmin) > Math.Abs(froot)) + { // Attempt inverse quadratic interpolation - s = froot / fxmin; - if (Close(xmin, xmax)) + double s = froot / fxmin; + double p; + double q; + if (xmin.AlmostEqual(xmax)) { p = 2.0 * xMid * s; q = 1.0 - s; @@ -67,22 +68,27 @@ namespace MathNet.Numerics.RootFinding else { q = fxmin / fxmax; - r = froot / fxmax; + double r = froot / fxmax; p = s * (2.0 * xMid * q * (q - r) - (root - xmin) * (r - 1.0)); q = (q - 1.0) * (r - 1.0) * (s - 1.0); } - if (p > 0.0) q = -q; // Check whether in bounds + + if (p > 0.0) + { + // Check whether in bounds + q = -q; + } p = Math.Abs(p); - min1 = 3.0 * xMid * q - Math.Abs(xAcc1 * q); - min2 = Math.Abs(e * q); - if (2.0 * p < Math.Min(min1, min2)) + if (2.0 * p < Math.Min(3.0 * xMid * q - Math.Abs(xAcc1 * q), Math.Abs(e * q))) { - e = d; // Accept interpolation + // Accept interpolation + e = d; d = p / q; } else { - d = xMid; // Interpolation failed, use bisection + // Interpolation failed, use bisection + d = xMid; e = d; } } @@ -95,9 +101,13 @@ namespace MathNet.Numerics.RootFinding xmin = root; fxmin = froot; if (Math.Abs(d) > xAcc1) + { root += d; + } else + { root += Sign(xAcc1, xMid); + } froot = f(root); } @@ -111,10 +121,5 @@ namespace MathNet.Numerics.RootFinding { return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); } - - static bool Close(double d1, double d2) - { - return Math.Abs(d1 - d2) <= double.Epsilon; - } } }