From 915c0b2c822d59d59099644b821a738a49d9aa42 Mon Sep 17 00:00:00 2001 From: rfellers Date: Tue, 1 Sep 2020 10:29:09 -0500 Subject: [PATCH] adds code from port, passing first test --- .../OptimizationTests/BrentMinimizerTests.cs | 2 +- .../GoldenSectionMinimizerTests.cs | 2 +- src/Numerics/Optimization/BrentMinimizer.cs | 128 +++++++++++++++++- 3 files changed, 128 insertions(+), 4 deletions(-) diff --git a/src/Numerics.Tests/OptimizationTests/BrentMinimizerTests.cs b/src/Numerics.Tests/OptimizationTests/BrentMinimizerTests.cs index e9decfed..98e43f34 100644 --- a/src/Numerics.Tests/OptimizationTests/BrentMinimizerTests.cs +++ b/src/Numerics.Tests/OptimizationTests/BrentMinimizerTests.cs @@ -52,7 +52,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests var algorithm = new BrentMinimizer(1e-5, 1000); var f1 = new Func(x => (x - 3) * (x - 3)); var obj = ObjectiveFunction.ScalarValue(f1); - var r1 = algorithm.FindMinimum(obj, -5, 5); + var r1 = algorithm.FindMinimum(obj, -2, 2); Assert.That(Math.Abs(r1.MinimizingPoint - 3.0), Is.LessThan(1e-4)); } diff --git a/src/Numerics.Tests/OptimizationTests/GoldenSectionMinimizerTests.cs b/src/Numerics.Tests/OptimizationTests/GoldenSectionMinimizerTests.cs index c4329beb..f81ec883 100644 --- a/src/Numerics.Tests/OptimizationTests/GoldenSectionMinimizerTests.cs +++ b/src/Numerics.Tests/OptimizationTests/GoldenSectionMinimizerTests.cs @@ -53,7 +53,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests var algorithm = new GoldenSectionMinimizer(1e-5, 1000); var f1 = new Func(x => (x - 3)*(x - 3)); var obj = ObjectiveFunction.ScalarValue(f1); - var r1 = algorithm.FindMinimum(obj, -5, 5); + var r1 = algorithm.FindMinimum(obj, -1, 1); Assert.That(Math.Abs(r1.MinimizingPoint - 3.0), Is.LessThan(1e-4)); } diff --git a/src/Numerics/Optimization/BrentMinimizer.cs b/src/Numerics/Optimization/BrentMinimizer.cs index 5b243fa3..792ce72b 100644 --- a/src/Numerics/Optimization/BrentMinimizer.cs +++ b/src/Numerics/Optimization/BrentMinimizer.cs @@ -27,6 +27,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using MathNet.Numerics.Optimization.ObjectiveFunctions; using System; namespace MathNet.Numerics.Optimization @@ -53,9 +54,132 @@ namespace MathNet.Numerics.Optimization return Minimum(objective, lowerBound, upperBound, XTolerance, MaximumIterations, MaximumExpansionSteps, LowerExpansionFactor, UpperExpansionFactor); } - public static ScalarMinimizationResult Minimum(IScalarObjectiveFunction objective, double lowerBound, double upperBound, double xTolerance = 1e-5, int maxIterations = 1000, int maxExpansionSteps = 10, double lowerExpansionFactor = 2.0, double upperExpansionFactor = 2.0) + public static ScalarMinimizationResult Minimum(IScalarObjectiveFunction objective, double lowerBound, double upperBound, double xTolerance = 1e-5, + int maxIterations = 1000, int maxExpansionSteps = 10, double lowerExpansionFactor = 2.0, double upperExpansionFactor = 2.0) { - return null; + int maxfun = maxIterations; + + if (lowerBound > upperBound) + throw new OptimizationException("Lower bound must be lower than upper bound."); + + double sqrt_eps = Math.Sqrt(2.2e-16); + + // This is not the golden_mean, but golden angle. Not sure why. + // https://en.wikipedia.org/wiki/Golden_angle + double golden_angle = 0.5 * (3.0 - Math.Sqrt(5.0)); + + double a = lowerBound; + double b = upperBound; + double fulc = a + golden_angle * (b - a); + + double nfc = fulc, xf = fulc; + double rat = 0.0, e = 0.0; + double x = xf; + var evaluation = objective.Evaluate(x); + double fx = evaluation.Value; + int num = 1; + + double fu = double.PositiveInfinity; + double ffulc = fx, fnfc = fx; + double xm = 0.5 * (a + b); + double tol1 = sqrt_eps * Math.Abs(xf) + xTolerance / 3.0; + double tol2 = 2.0 * tol1; + + while (Math.Abs(xf - xm) > (tol2 - 0.5 * (b - a))) + { + bool golden = true; + + // Check for parabolic fit + if (Math.Abs(e) > tol1) + { + golden = false; + double r = (xf - nfc) * (fx - ffulc); + double q = (xf - fulc) * (fx - fnfc); + double p = (xf - fulc) * q - (xf - nfc) * r; + q = 2.0 * (q - r); + if (q > 0.0) + p = -p; + q = Math.Abs(q); + r = e; + e = rat; + + // Check for acceptability of parabola + if ((Math.Abs(p) < Math.Abs(0.5 * q * r)) && (p > q * (a - xf)) && (p < q * (b - xf))) + { + rat = (p + 0.0) / q; + x = xf + rat; + + if (((x - a) < tol2) || ((b - x) < tol2)) + { + int si_2 = Math.Sign(xm - xf) + ((xm - xf) == 0 ? 1 : 0); + rat = tol1 * si_2; + } + } + else // do a golden-section step + golden = true; + } + + if (golden) // do a golden-section step + { + if (xf >= xm) + e = a - xf; + else + e = b - xf; + rat = golden_angle * e; + } + + int si = Math.Sign(rat) + (rat == 0 ? 1 : 0); + x = xf + si * Math.Max(Math.Abs(rat), tol1); + + evaluation = objective.Evaluate(x); + fu = evaluation.Value; + num += 1; + + if (fu <= fx) + { + if (x >= xf) + a = xf; + else + b = xf; + + fulc = nfc; ffulc = fnfc; + nfc = xf; fnfc = fx; + xf = x; fx = fu; + } + else + { + if (x < xf) + a = x; + else + b = x; + + if ((fu <= fnfc) || (nfc == xf)) + { + fulc = nfc; ffulc = fnfc; + nfc = x; fnfc = fu; + } + else if ((fu <= ffulc) || (fulc == xf) || (fulc == nfc)) + { + fulc = x; ffulc = fu; + } + } + + xm = 0.5 * (a + b); + tol1 = sqrt_eps * Math.Abs(xf) + xTolerance / 3.0; + tol2 = 2.0 * tol1; + + if (num >= maxfun) + break; + } + + var exitCondition = ExitCondition.BoundTolerance; + + if (num >= maxfun) + exitCondition = ExitCondition.ExceedIterations; + else if (double.IsNaN(xf) || double.IsNaN(fx) || double.IsNaN(fu)) + exitCondition = ExitCondition.InvalidValues; + + return new ScalarMinimizationResult(new ScalarValueObjectiveFunctionEvaluation(xf, fx), num, exitCondition); } static void ValueChecker(double value, double point)