From f621a940d12a2cc41cf75f990fb6259202f8a0b0 Mon Sep 17 00:00:00 2001 From: Scott Stephens Date: Fri, 23 Sep 2016 16:34:20 -0500 Subject: [PATCH] Optimization: revamp testing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Changed file and class names from TestClassName to ClassNameTests to better conform with other parts of Math.Net * Added some organizational abstractions to simplify adding test cases. * Added test cases based on several functions from Testing Unconstrained Optimization Software Jorge J. Moré, Burton S. Garbow, Kenneth E. Hillstrom ACM Transactions on Mathematical Software, Vol 7, No. 1, March 1981, Pages 17-41. * At this point, most test functions are low-dimensional, well-scaled, unimodal functions. * There are test failures for Conjugate Gradient and Newton unconstrained minimizers. Currently I believe these to be failures of the algorithm, rather than failures of implementation, but it hasn't been investigated thoroughly. --- ...gsBMinimizer.cs => BfgsBMinimizerTests.cs} | 55 +++- ...BfgsMinimizer.cs => BfgsMinimizerTests.cs} | 55 +++- .../ConjugateGradientMinimizerTests.cs | 85 +++++ ...izer.cs => GoldenSectionMinimizerTests.cs} | 2 +- ...onMinimizer.cs => NewtonMinimizerTests.cs} | 54 ++- ...Function.cs => RosenbrockFunctionTests.cs} | 2 +- .../TestConjugateGradientMinimizer.cs | 33 -- .../OptimizationTests/TestFunctionAdapters.cs | 54 +++ .../OptimizationTests/TestFunctionTests.cs | 309 ++++++++++++++++++ .../TestFunctions/BaseTestFunction.cs | 125 +++++++ .../TestFunctions/BealeFunction.cs | 97 ++++++ .../TestFunctions/BrownAndDennisFunction.cs | 114 +++++++ .../TestFunctions/BrownBadlyScaledFunction.cs | 131 ++++++++ .../FreudensteinAndRothFunction.cs | 98 ++++++ .../TestFunctions/HelicalValleyFunction.cs | 178 ++++++++++ .../TestFunctions/ITestFunction.cs | 69 ++++ .../JennrichAndSampsonFunction.cs | 102 ++++++ .../TestFunctions/MeyerFunction.cs | 99 ++++++ .../PowellBadlyScaledFunction.cs | 111 +++++++ .../TestFunctions/PowellSingularFunction.cs | 141 ++++++++ .../TestFunctions/RosenbrockFunction2.cs | 156 +++++++++ .../TestFunctions/WoodFunction.cs | 164 ++++++++++ src/UnitTests/UnitTests.csproj | 27 +- 23 files changed, 2217 insertions(+), 44 deletions(-) rename src/UnitTests/OptimizationTests/{TestBfgsBMinimizer.cs => BfgsBMinimizerTests.cs} (75%) rename src/UnitTests/OptimizationTests/{TestBfgsMinimizer.cs => BfgsMinimizerTests.cs} (62%) create mode 100644 src/UnitTests/OptimizationTests/ConjugateGradientMinimizerTests.cs rename src/UnitTests/OptimizationTests/{TestGoldenSectionMinimizer.cs => GoldenSectionMinimizerTests.cs} (95%) rename src/UnitTests/OptimizationTests/{TestNewtonMinimizer.cs => NewtonMinimizerTests.cs} (68%) rename src/UnitTests/OptimizationTests/{TestRosenbrockFunction.cs => RosenbrockFunctionTests.cs} (98%) delete mode 100644 src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctionAdapters.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctionTests.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/BaseTestFunction.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/BealeFunction.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/BrownAndDennisFunction.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/BrownBadlyScaledFunction.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/FreudensteinAndRothFunction.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/HelicalValleyFunction.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/ITestFunction.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/JennrichAndSampsonFunction.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/MeyerFunction.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/PowellBadlyScaledFunction.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/PowellSingularFunction.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/RosenbrockFunction2.cs create mode 100644 src/UnitTests/OptimizationTests/TestFunctions/WoodFunction.cs diff --git a/src/UnitTests/OptimizationTests/TestBfgsBMinimizer.cs b/src/UnitTests/OptimizationTests/BfgsBMinimizerTests.cs similarity index 75% rename from src/UnitTests/OptimizationTests/TestBfgsBMinimizer.cs rename to src/UnitTests/OptimizationTests/BfgsBMinimizerTests.cs index da70b9d1..785ddb81 100644 --- a/src/UnitTests/OptimizationTests/TestBfgsBMinimizer.cs +++ b/src/UnitTests/OptimizationTests/BfgsBMinimizerTests.cs @@ -32,11 +32,16 @@ using System; using MathNet.Numerics.LinearAlgebra.Double; using MathNet.Numerics.Optimization; using NUnit.Framework; +using System.Linq; +using System.Text; +using System.Collections.Generic; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using System.Collections; namespace MathNet.Numerics.UnitTests.OptimizationTests { [TestFixture] - public class TestBfgsBMinimizer + public class BfgsBMinimizerTests { [Test] public void FindMinimum_Rosenbrock_Easy() @@ -149,6 +154,54 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests Assert.LessOrEqual(result.MinimizingPoint[0],upperBound[0]); Assert.LessOrEqual(result.MinimizingPoint[1],upperBound[1]); } + + [Test] + [TestCaseSource(typeof(MghTestCaseEnumerator))] + public void Mgh_Tests(TestFunctions.TestCase test_case) + { + var obj = new MghObjectiveFunction(test_case.Function, true, true); + var solver = new BfgsBMinimizer(1e-8, 1e-8, 1e-8, 1000); + + var result = solver.FindMinimum(obj, test_case.LowerBound, test_case.UpperBound, test_case.InitialGuess); + + if (test_case.MinimizingPoint != null) + { + Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3)); + } + + var val1 = result.FunctionInfoAtMinimum.Value; + var val2 = test_case.MinimalValue; + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + var abs_err = Math.Abs(val1 - val2); + var rel_err = abs_err / abs_min; + var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3); + Assert.That(success, "Minimal function value is not as expected."); + } + + private class MghTestCaseEnumerator : IEnumerable + { + public IEnumerator GetEnumerator() + { + return + RosenbrockFunction2.TestCases + .Concat(BealeFunction.TestCases) + .Concat(HelicalValleyFunction.TestCases) + .Concat(MeyerFunction.TestCases) + .Concat(PowellSingularFunction.TestCases) + .Concat(WoodFunction.TestCases) + .Concat(BrownAndDennisFunction.TestCases) + .Where(x => x.IsBounded) + .Select(x => new TestCaseData(x) + .SetName(x.FullName) + ) + .GetEnumerator(); + } + + IEnumerator IEnumerable.GetEnumerator() + { + return this.GetEnumerator(); + } + } } } diff --git a/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs b/src/UnitTests/OptimizationTests/BfgsMinimizerTests.cs similarity index 62% rename from src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs rename to src/UnitTests/OptimizationTests/BfgsMinimizerTests.cs index b2300d35..be4e4b97 100644 --- a/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs +++ b/src/UnitTests/OptimizationTests/BfgsMinimizerTests.cs @@ -1,12 +1,17 @@ using System; +using System.Linq; using MathNet.Numerics.LinearAlgebra.Double; using MathNet.Numerics.Optimization; using NUnit.Framework; +using System.Text; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using System.Collections.Generic; +using System.Collections; namespace MathNet.Numerics.UnitTests.OptimizationTests { [TestFixture] - public class TestBfgsMinimizer + public class BfgsMinimizerTests { [Test] public void FindMinimum_Rosenbrock_Easy() @@ -73,5 +78,53 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests Assert.That(Math.Abs(result.MinimizingPoint[0] - BigRosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); Assert.That(Math.Abs(result.MinimizingPoint[1] - BigRosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); } + + private class MghTestCaseEnumerator : IEnumerable + { + public IEnumerator GetEnumerator() + { + return + RosenbrockFunction2.TestCases + .Concat(BealeFunction.TestCases) + .Concat(HelicalValleyFunction.TestCases) + .Concat(MeyerFunction.TestCases) + .Concat(PowellSingularFunction.TestCases) + .Concat(WoodFunction.TestCases) + .Concat(BrownAndDennisFunction.TestCases) + .Where(x => x.IsUnbounded) + .Select(x => new TestCaseData(x) + .SetName(x.FullName) + ) + .GetEnumerator(); + } + + IEnumerator IEnumerable.GetEnumerator() + { + return this.GetEnumerator(); + } + } + + [Test] + [TestCaseSource(typeof(MghTestCaseEnumerator))] + public void Mgh_Tests(TestFunctions.TestCase test_case) + { + var obj = new MghObjectiveFunction(test_case.Function, true, true); + var solver = new BfgsMinimizer(1e-8, 1e-8, 1e-8, 1000); + + var result = solver.FindMinimum(obj, test_case.InitialGuess); + + if (test_case.MinimizingPoint != null) + { + Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3)); + } + + var val1 = result.FunctionInfoAtMinimum.Value; + var val2 = test_case.MinimalValue; + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + var abs_err = Math.Abs(val1 - val2); + var rel_err = abs_err / abs_min; + var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3); + Assert.That(success, "Minimal function value is not as expected."); + } } } diff --git a/src/UnitTests/OptimizationTests/ConjugateGradientMinimizerTests.cs b/src/UnitTests/OptimizationTests/ConjugateGradientMinimizerTests.cs new file mode 100644 index 00000000..4a13ba65 --- /dev/null +++ b/src/UnitTests/OptimizationTests/ConjugateGradientMinimizerTests.cs @@ -0,0 +1,85 @@ +using System; +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Optimization; +using NUnit.Framework; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using System.Collections; +using System.Collections.Generic; +using System.Linq; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + public class ConjugateGradientMinimizerTests + { + [Test] + public void FindMinimum_Rosenbrock_Easy() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new ConjugateGradientMinimizer(1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[]{1.2,1.2})); + + Assert.That(Math.Abs(result.MinimizingPoint[0]-1.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Hard() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new ConjugateGradientMinimizer(1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2, 1.0 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + private class MghTestCaseEnumerator : IEnumerable + { + public IEnumerator GetEnumerator() + { + return + RosenbrockFunction2.TestCases + .Concat(BealeFunction.TestCases) + .Concat(HelicalValleyFunction.TestCases) + .Concat(MeyerFunction.TestCases) + .Concat(PowellSingularFunction.TestCases) + .Concat(WoodFunction.TestCases) + .Concat(BrownAndDennisFunction.TestCases) + .Where(x => x.IsUnbounded) + .Select(x => new TestCaseData(x) + .SetName(x.FullName) + ) + .GetEnumerator(); + } + + IEnumerator IEnumerable.GetEnumerator() + { + return this.GetEnumerator(); + } + } + + [Test] + [TestCaseSource(typeof(MghTestCaseEnumerator))] + public void Mgh_Tests(TestFunctions.TestCase test_case) + { + var obj = new MghObjectiveFunction(test_case.Function, true, true); + var solver = new ConjugateGradientMinimizer(1e-8, 1000); + + var result = solver.FindMinimum(obj, test_case.InitialGuess); + + if (test_case.MinimizingPoint != null) + { + Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3)); + } + + var val1 = result.FunctionInfoAtMinimum.Value; + var val2 = test_case.MinimalValue; + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + var abs_err = Math.Abs(val1 - val2); + var rel_err = abs_err / abs_min; + var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3); + Assert.That(success, "Minimal function value is not as expected."); + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestGoldenSectionMinimizer.cs b/src/UnitTests/OptimizationTests/GoldenSectionMinimizerTests.cs similarity index 95% rename from src/UnitTests/OptimizationTests/TestGoldenSectionMinimizer.cs rename to src/UnitTests/OptimizationTests/GoldenSectionMinimizerTests.cs index a85fa7f1..e04088d9 100644 --- a/src/UnitTests/OptimizationTests/TestGoldenSectionMinimizer.cs +++ b/src/UnitTests/OptimizationTests/GoldenSectionMinimizerTests.cs @@ -5,7 +5,7 @@ using NUnit.Framework; namespace MathNet.Numerics.UnitTests.OptimizationTests { [TestFixture] - public class TestGoldenSectionMinimizer + public class GoldenSectionMinimizerTests { [Test] public void Test_Works() diff --git a/src/UnitTests/OptimizationTests/TestNewtonMinimizer.cs b/src/UnitTests/OptimizationTests/NewtonMinimizerTests.cs similarity index 68% rename from src/UnitTests/OptimizationTests/TestNewtonMinimizer.cs rename to src/UnitTests/OptimizationTests/NewtonMinimizerTests.cs index 772c07ad..efccfa14 100644 --- a/src/UnitTests/OptimizationTests/TestNewtonMinimizer.cs +++ b/src/UnitTests/OptimizationTests/NewtonMinimizerTests.cs @@ -3,6 +3,10 @@ using MathNet.Numerics.LinearAlgebra.Double; using MathNet.Numerics.Optimization; using MathNet.Numerics.Optimization.ObjectiveFunctions; using NUnit.Framework; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using System.Collections.Generic; +using System.Collections; +using System.Linq; namespace MathNet.Numerics.UnitTests.OptimizationTests { @@ -51,7 +55,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests } [TestFixture] - public class TestNewtonMinimizer + public class NewtonMinimizerTests { [Test] public void FindMinimum_Rosenbrock_Easy() @@ -118,5 +122,53 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3)); Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); } + + private class MghTestCaseEnumerator : IEnumerable + { + public IEnumerator GetEnumerator() + { + return + RosenbrockFunction2.TestCases + .Concat(BealeFunction.TestCases) + .Concat(HelicalValleyFunction.TestCases) + .Concat(MeyerFunction.TestCases) + .Concat(PowellSingularFunction.TestCases) + .Concat(WoodFunction.TestCases) + .Concat(BrownAndDennisFunction.TestCases) + .Where(x => x.IsUnbounded) + .Select(x => new TestCaseData(x) + .SetName(x.FullName) + ) + .GetEnumerator(); + } + + IEnumerator IEnumerable.GetEnumerator() + { + return this.GetEnumerator(); + } + } + + [Test] + [TestCaseSource(typeof(MghTestCaseEnumerator))] + public void Mgh_Tests(TestFunctions.TestCase test_case) + { + var obj = new MghObjectiveFunction(test_case.Function, true, true); + var solver = new NewtonMinimizer(1e-8, 1000, useLineSearch: false); + + var result = solver.FindMinimum(obj, test_case.InitialGuess); + + if (test_case.MinimizingPoint != null) + { + Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3)); + } + + var val1 = result.FunctionInfoAtMinimum.Value; + var val2 = test_case.MinimalValue; + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + var abs_err = Math.Abs(val1 - val2); + var rel_err = abs_err / abs_min; + var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3); + Assert.That(success, "Minimal function value is not as expected."); + } } } diff --git a/src/UnitTests/OptimizationTests/TestRosenbrockFunction.cs b/src/UnitTests/OptimizationTests/RosenbrockFunctionTests.cs similarity index 98% rename from src/UnitTests/OptimizationTests/TestRosenbrockFunction.cs rename to src/UnitTests/OptimizationTests/RosenbrockFunctionTests.cs index 14a0f9f7..b1b0e5ae 100644 --- a/src/UnitTests/OptimizationTests/TestRosenbrockFunction.cs +++ b/src/UnitTests/OptimizationTests/RosenbrockFunctionTests.cs @@ -5,7 +5,7 @@ using NUnit.Framework; namespace MathNet.Numerics.UnitTests.OptimizationTests { [TestFixture] - class TestRosenbrockFunction + class RosenbrockFunctionTests { [Test] public void TestGradient() diff --git a/src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs b/src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs deleted file mode 100644 index bd1927e5..00000000 --- a/src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs +++ /dev/null @@ -1,33 +0,0 @@ -using System; -using MathNet.Numerics.LinearAlgebra.Double; -using MathNet.Numerics.Optimization; -using NUnit.Framework; - -namespace MathNet.Numerics.UnitTests.OptimizationTests -{ - [TestFixture] - public class TestConjugateGradientMinimizer - { - [Test] - public void FindMinimum_Rosenbrock_Easy() - { - var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); - var solver = new ConjugateGradientMinimizer(1e-5, 1000); - var result = solver.FindMinimum(obj, new DenseVector(new[]{1.2,1.2})); - - Assert.That(Math.Abs(result.MinimizingPoint[0]-1.0), Is.LessThan(1e-3)); - Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); - } - - [Test] - public void FindMinimum_Rosenbrock_Hard() - { - var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); - var solver = new ConjugateGradientMinimizer(1e-5, 1000); - var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2, 1.0 })); - - Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3)); - Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); - } - } -} diff --git a/src/UnitTests/OptimizationTests/TestFunctionAdapters.cs b/src/UnitTests/OptimizationTests/TestFunctionAdapters.cs new file mode 100644 index 00000000..8f20f954 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctionAdapters.cs @@ -0,0 +1,54 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.Optimization; +using MathNet.Numerics.Optimization.ObjectiveFunctions; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using MathNet.Numerics.LinearAlgebra.Double; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + public class MghObjectiveFunction : LazyObjectiveFunctionBase + { + private ITestFunction TestFunction; + + public MghObjectiveFunction(ITestFunction testFunction, bool use_gradient, bool use_hessian) + : base(use_gradient, use_hessian) + { + this.TestFunction = testFunction; + } + + public override IObjectiveFunction CreateNew() + { + return new MghObjectiveFunction(this.TestFunction, this.IsGradientSupported, this.IsHessianSupported); + } + + protected override void EvaluateValue() + { + this.Value = this.TestFunction.SsqValue(this.Point); + } + + protected override void EvaluateGradient() + { + if (this.IsGradientSupported) + { + if (this._gradientValue == null) + this.Gradient = new DenseVector(this.TestFunction.ParameterDimension); + this.TestFunction.SsqGradientByRef(this.Point, _gradientValue); + } + } + + protected override void EvaluateHessian() + { + if (this.IsHessianSupported) + { + if (this._hessianValue == null) + this.Hessian = new DenseMatrix(this.TestFunction.ParameterDimension, this.TestFunction.ParameterDimension); + this.TestFunction.SsqHessianByRef(this.Point, _hessianValue); + } + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctionTests.cs b/src/UnitTests/OptimizationTests/TestFunctionTests.cs new file mode 100644 index 00000000..e3108615 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctionTests.cs @@ -0,0 +1,309 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using NUnit.Framework; +using MathNet.Numerics.LinearAlgebra; +using System.Collections; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + public class TestFunctionTests + { + private static IEnumerable MghCases + { + get + { + return Enumerable.Empty() + .Concat(RosenbrockFunction2.TestCases) + .Concat(BealeFunction.TestCases) + .Concat(HelicalValleyFunction.TestCases) + .Concat(MeyerFunction.TestCases) + .Concat(PowellSingularFunction.TestCases) + .Concat(WoodFunction.TestCases) + .Concat(BrownAndDennisFunction.TestCases); + } + } + + private class MghCaseEnumerator : IEnumerable + { + public string CategoryName { get; protected set; } + + public MghCaseEnumerator(string category_name) + { + this.CategoryName = category_name; + } + + public virtual IEnumerator GetEnumerator() + { + return MghCases + .Select(x => + new TestCaseData(x) + .SetName($"{x.FullName} {this.CategoryName}") + ).GetEnumerator(); + } + + IEnumerator IEnumerable.GetEnumerator() + { + return this.GetEnumerator(); + } + } + + [Test] + public void Smoke_Construction() + { + var c = new TestCase() + { + InitialGuess = new double[] { 1, 2, 3 }, + MinimizingPoint = new double[] { 1, 1, 1 }, + MinimalValue = 0 + }; + } + + private class ValueAtMinimumSource : MghCaseEnumerator + { + public ValueAtMinimumSource() : base("ValueAtMinimum") { } + + public override IEnumerator GetEnumerator() + { + return MghCases + .Where(x => x.MinimizingPoint != null) + .Select(x => + new TestCaseData(x) + .SetName($"{x.FullName} {this.CategoryName}") + ) + .GetEnumerator(); + } + } + + [Test] + [TestCaseSource(typeof(ValueAtMinimumSource))] + public void ValueAtMinimum(TestFunctions.TestCase test_case) + { + if (test_case.MinimizingPoint != null) + { + var value_at_minimum = test_case.Function.SsqValue(test_case.MinimizingPoint); + Assert.That( + Math.Abs(value_at_minimum - test_case.MinimalValue) < 1e-3, + $"Function value at minimum not as expected." + ); + } + } + + private class GradientAtStartSource : MghCaseEnumerator + { + public GradientAtStartSource() : base("GradientAtStart") { } + } + + [Test] + [TestCaseSource(typeof(GradientAtStartSource))] + public void GradientAtStart(TestFunctions.TestCase test_case) + { + var a_grad = test_case.Function.SsqGradient(test_case.InitialGuess); + var fd_grad = Vector.Build.Dense(test_case.Function.ParameterDimension, 0.0); + + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + var h = 1e-6; + + var bump_up = test_case.InitialGuess.Clone(); + bump_up[ii] += h; + var bump_down = test_case.InitialGuess.Clone(); + bump_down[ii] -= h; + + var up_val = test_case.Function.SsqValue(bump_up); + var down_val = test_case.Function.SsqValue(bump_down); + + fd_grad[ii] = 0.5 * (up_val - down_val) / h; + } + + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + var val1 = a_grad[ii]; + var val2 = fd_grad[ii]; + var min_abs_val = Math.Min(Math.Abs(val1), Math.Abs(val2)); + if (min_abs_val <= 1) + Assert.That(Math.Abs(val1 - val2) < 1e-3, $"Problem with gradient value at start point."); + else + Assert.That(Math.Abs(val1 - val2) / min_abs_val < 1e-3, $"Problem with gradient value at start point."); + } + } + + private class HessianAtStartSource : MghCaseEnumerator + { + public HessianAtStartSource() : base("HessianAtStart") { } + + public override IEnumerator GetEnumerator() + { + return MghCases + .Where(x => x.MinimizingPoint != null) + .Select(x => + new TestCaseData(x) + .SetName($"{x.FullName} {this.CategoryName}") + ) + .GetEnumerator(); + } + } + [Test] + [TestCaseSource(typeof(HessianAtStartSource))] + public void HessianAtStart(TestFunctions.TestCase test_case) + { + var a_hess = test_case.Function.SsqHessian(test_case.InitialGuess); + var fd_hess = Matrix.Build.Dense(test_case.Function.ParameterDimension, test_case.Function.ParameterDimension); + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + for (int jj = 0; jj < test_case.Function.ParameterDimension; ++jj) + { + var h1 = 1e-3 * Math.Max(1.0, Math.Abs(test_case.InitialGuess[ii])); + var h2 = 1e-3 * Math.Max(1.0, Math.Abs(test_case.InitialGuess[jj])); + + var bump_uu = test_case.InitialGuess.Clone(); + bump_uu[ii] += h1; + bump_uu[jj] += h2; + + var bump_dd = test_case.InitialGuess.Clone(); + bump_dd[ii] -= h1; + bump_dd[jj] -= h2; + + var bump_ud = test_case.InitialGuess.Clone(); + bump_ud[ii] += h1; + bump_ud[jj] -= h2; + + var bump_du = test_case.InitialGuess.Clone(); + bump_du[ii] -= h1; + bump_du[jj] += h2; + + var val_uu = test_case.Function.SsqValue(bump_uu); + var val_dd = test_case.Function.SsqValue(bump_dd); + var val_ud = test_case.Function.SsqValue(bump_ud); + var val_du = test_case.Function.SsqValue(bump_du); + + fd_hess[ii, jj] = (val_uu - val_ud + val_dd - val_du) / (4 * h1 * h2); + } + } + + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + for (int jj = 0; jj < test_case.Function.ParameterDimension; ++jj) + { + var val1 = fd_hess[ii, jj]; + var val2 = a_hess[ii, jj]; + + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + if (abs_min <= 1) + { + Assert.That(Math.Abs(val1 - val2) < 1e-3, $"Problem with hessian at start point."); + } + else + { + Assert.That(Math.Abs(val1 - val2) / abs_min < 0.05, $"Problem with hessian at start point."); + } + } + } + } + + private class ItemGradientAtStartSource : MghCaseEnumerator + { + public ItemGradientAtStartSource() : base("ItemGradientAtStart") { } + } + + [Test] + [TestCaseSource(typeof(ItemGradientAtStartSource))] + public void ItemGradientAtStart(TestFunctions.TestCase test_case) + { + for (var item_index = 0; item_index < test_case.Function.ItemDimension; ++item_index) + { + + var a_grad = test_case.Function.ItemGradient(test_case.InitialGuess, item_index); + var h = 1e-4; + var fd_grad = Vector.Build.Dense(test_case.Function.ParameterDimension, 0.0); + + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + var bump_up = test_case.InitialGuess.Clone(); + bump_up[ii] += h; + var bump_down = test_case.InitialGuess.Clone(); + bump_down[ii] -= h; + + var up_val = test_case.Function.ItemValue(bump_up, item_index); + var down_val = test_case.Function.ItemValue(bump_down, item_index); + + fd_grad[ii] = 0.5 * (up_val - down_val) / h; + } + + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + Assert.That(Math.Abs(fd_grad[ii] - a_grad[ii]) < 1e-3, $"Failed for parameter {ii}"); + } + } + } + + private class ItemHessianAtStartSource : MghCaseEnumerator + { + public ItemHessianAtStartSource() : base("ItemHessianAtStart") { } + } + + [Test] + [TestCaseSource(typeof(ItemHessianAtStartSource))] + public void ItemHessianAtStart(TestFunctions.TestCase test_case) + { + for (var item_index = 0; item_index < test_case.Function.ItemDimension; ++item_index) + { + var a_hess = test_case.Function.ItemHessian(test_case.InitialGuess, item_index); + var h = 1e-4; + var fd_hess = Matrix.Build.Dense(test_case.Function.ParameterDimension, test_case.Function.ParameterDimension); + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + for (int jj = 0; jj < test_case.Function.ParameterDimension; ++jj) + { + var bump_uu = test_case.InitialGuess.Clone(); + bump_uu[ii] += h; + bump_uu[jj] += h; + + var bump_dd = test_case.InitialGuess.Clone(); + bump_dd[ii] -= h; + bump_dd[jj] -= h; + + var bump_ud = test_case.InitialGuess.Clone(); + bump_ud[ii] += h; + bump_ud[jj] -= h; + + var bump_du = test_case.InitialGuess.Clone(); + bump_du[ii] -= h; + bump_du[jj] += h; + + var val_uu = test_case.Function.ItemValue(bump_uu, item_index); + var val_dd = test_case.Function.ItemValue(bump_dd, item_index); + var val_ud = test_case.Function.ItemValue(bump_ud, item_index); + var val_du = test_case.Function.ItemValue(bump_du, item_index); + + fd_hess[ii, jj] = (val_uu - val_ud + val_dd - val_du) / (4 * h * h); + } + } + + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + for (int jj = 0; jj < test_case.Function.ParameterDimension; ++jj) + { + var val1 = fd_hess[ii, jj]; + var val2 = a_hess[ii, jj]; + + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + if (abs_min <= 1) + { + Assert.That(Math.Abs(val1 - val2) < 1e-3, $"Problem with hessian at start point."); + } + else + { + Assert.That(Math.Abs(val1 - val2) / abs_min < 0.05, $"Problem with hessian at start point."); + } + } + } + } + } + + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/BaseTestFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/BaseTestFunction.cs new file mode 100644 index 00000000..60da6c61 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/BaseTestFunction.cs @@ -0,0 +1,125 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public abstract class BaseTestFunction : ITestFunction + { + public abstract string Description { get; } + public abstract int ParameterDimension { get; } + public abstract int ItemDimension { get; } + + public abstract double ItemValue(Vector x, int itemIndex); + public abstract void ItemGradientByRef(Vector x, int itemIndex, Vector output); + public abstract void ItemHessianByRef(Vector x, int itemIndex, Matrix output); + + public virtual Vector ItemGradient(Vector x, int itemIndex) + { + var output = new LinearAlgebra.Double.DenseVector(this.ParameterDimension); + this.ItemGradientByRef(x, itemIndex, output); + return output; + } + + public virtual Matrix ItemHessian(Vector x, int itemIndex) + { + var output = new LinearAlgebra.Double.DenseMatrix(this.ParameterDimension, this.ParameterDimension); + this.ItemHessianByRef(x, itemIndex, output); + return output; + } + + + public virtual void JacobianbyRef(Vector x, Matrix output) + { + for (int ii = 0; ii < this.ItemDimension; ++ii) + { + var grad = this.ItemGradient(x, ii); + output.SetRow(ii, grad); + } + } + + public virtual Matrix Jacobian(Vector x) + { + var output = new LinearAlgebra.Double.DenseMatrix(this.ItemDimension, this.ParameterDimension); + this.JacobianbyRef(x, output); + return output; + } + + public virtual void SsqGradientByRef(Vector x, Vector output) + { + if (output.Count != this.ParameterDimension) + throw new ArgumentException($"Output vector must match parameter dimension of function; expected {this.ParameterDimension}, got {output.Count}."); + + for (int jj = 0; jj < this.ParameterDimension; ++jj) + output[jj] = 0.0; + var tmp_grad = new LinearAlgebra.Double.DenseVector(this.ParameterDimension); + double tmp_value = 0.0; + + for (int ii = 0; ii < this.ItemDimension; ++ii) + { + tmp_value = this.ItemValue(x, ii); + this.ItemGradientByRef(x, ii, tmp_grad); + for (int jj = 0; jj < this.ParameterDimension; ++jj) + output[jj] += 2 * tmp_value * tmp_grad[jj]; + } + } + + public virtual Vector SsqGradient(Vector x) + { + var output = new LinearAlgebra.Double.DenseVector(this.ParameterDimension); + this.SsqGradientByRef(x, output); + return output; + } + + public virtual void SsqHessianByRef(Vector x, Matrix output) + { + if (output.RowCount != this.ParameterDimension || output.ColumnCount != this.ParameterDimension) + throw new ArgumentException($"Output matrix must match parameter dimension of function; expected {this.ParameterDimension}x{this.ParameterDimension}, got {output.RowCount}x{output.ColumnCount}."); + + for (int ii = 0; ii < this.ParameterDimension; ++ii) + for (int jj = 0; jj < this.ParameterDimension; ++jj) + output[ii,jj] = 0.0; + + var tmp_grad = new LinearAlgebra.Double.DenseVector(this.ParameterDimension); + var tmp_hess = new LinearAlgebra.Double.DenseMatrix(this.ParameterDimension, this.ParameterDimension); + double tmp_value = 0.0; + + for (int ii = 0; ii < this.ItemDimension; ++ii) + { + tmp_value = this.ItemValue(x, ii); + this.ItemGradientByRef(x, ii, tmp_grad); + this.ItemHessianByRef(x, ii, tmp_hess); + for (int jj = 0; jj < this.ParameterDimension; ++jj) + { + for (int kk = 0; kk < this.ParameterDimension; ++kk) + { + var increment = 2 * (tmp_value * tmp_hess[jj, kk] + tmp_grad[jj] * tmp_grad[kk]); + output[jj, kk] += increment; + } + } + } + } + + public virtual Matrix SsqHessian(Vector x) + { + var output = new LinearAlgebra.Double.DenseMatrix(this.ParameterDimension, this.ParameterDimension); + this.SsqHessianByRef(x, output); + return output; + } + + public virtual double SsqValue(Vector x) + { + double ssq = 0.0; + for (int ii = 0; ii < this.ItemDimension; ++ii) + { + var tmp = this.ItemValue(x, ii); + ssq += tmp * tmp; + } + return ssq; + } + + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/BealeFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/BealeFunction.cs new file mode 100644 index 00000000..b55208dd --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/BealeFunction.cs @@ -0,0 +1,97 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class BealeFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new BealeFunction(), + InitialGuess = new double[] { 1, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 3, 0.5 }, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new BealeFunction(), + InitialGuess = new double[] { 1, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 3, 0.5 }, + LowerBound = new double[] { -1000, -1000}, + UpperBound = new double[] { 1000, 1000}, + CaseName = "loose bounds" + }; + yield return new TestCase() + { + Function = new BealeFunction(), + InitialGuess = new double[] { 1, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 3, 0.5 }, + LowerBound = new double[] { 0.6, 0.5 }, + UpperBound = new double[] { 10, 100 }, + CaseName = "tight bounds" + }; + } + } + + public BealeFunction() { } + + public override string Description + { + get + { + return "Beale fun (MGH #5)"; + } + } + + public override int ItemDimension + { + get + { + return 3; + } + } + + public override int ParameterDimension + { + get + { + return 2; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + int ii = itemIndex + 1; + output[0] = -1 + Math.Pow(x[1], ii); + output[1] = ii * x[0] * Math.Pow(x[1], ii - 1); + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + int ii = itemIndex + 1; + output[0, 0] = 0; + output[0, 1] = ii * Math.Pow(x[1], ii - 1); + output[1, 0] = ii * Math.Pow(x[1], ii - 1); + output[1, 1] = (ii - 1) * ii * x[0] * Math.Pow(x[1], ii - 2); + } + + private static readonly double[] y = { 1.5, 2.25, 2.625}; + + public override double ItemValue(Vector x, int itemIndex) + { + int ii = itemIndex + 1; + return y[itemIndex] - x[0] * (1 - Math.Pow(x[1], ii)); + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/BrownAndDennisFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/BrownAndDennisFunction.cs new file mode 100644 index 00000000..9e5b9c41 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/BrownAndDennisFunction.cs @@ -0,0 +1,114 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class BrownAndDennisFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new BrownAndDennisFunction(20), + InitialGuess = new double[] { 25, 5, -5, -1 }, + MinimalValue = 85822.2, + MinimizingPoint = null, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new BrownAndDennisFunction(20), + InitialGuess = new double[] { 25, 5, -5, -1 }, + MinimalValue = 85822.2, + MinimizingPoint = null, + LowerBound = new double[] { -1000, -1000, -1000, -1000 }, + UpperBound = new double[] {1000, 1000, 1000, 1000 }, + CaseName = "loose bounds" + }; + yield return new TestCase() + { + Function = new BrownAndDennisFunction(20), + InitialGuess = new double[] { 25, 5, -5, -1 }, + MinimalValue = 0.88860479e5, + MinimizingPoint = null, + LowerBound = new double[] { -10, 0, -100, -20 }, + UpperBound = new double[] { 100, 15, 0, 0.2 }, + CaseName = "tight bounds" + }; + } + } + + private readonly int _items; + + public BrownAndDennisFunction(int items) + { + if (items < 4) + throw new ArgumentException("items must be >= 4"); + _items = items; + } + + public override string Description + { + get + { + return "Brown & Dennis fun (MGH #16)"; + } + } + + public override int ItemDimension + { + get + { + return _items; + } + } + + public override int ParameterDimension + { + get + { + return 4; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + var ii = itemIndex + 1; + var t = ii / 5.0; + output[0] = 2 * (x[0] + t * x[1] - Math.Exp(t)); + output[1] = (2*ii/25.0) * (5 * x[0] + ii * x[1] - 5 * Math.Exp(t)); + output[2] = 2 * (x[2] + x[3] * Math.Sin(t) - Math.Cos(t)); + output[3] = 2 * Math.Sin(t) * (x[2] + Math.Sin(t) * x[3] - Math.Cos(t)); + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + for (int ii = 0; ii < 4; ++ii) + for (int jj = 0; jj < 4; ++jj) + output[ii, jj] = 0; + var i = itemIndex + 1; + var t = i / 5.0; + output[0, 0] = 2; + output[0, 1] = 2 * t; + output[1, 0] = 2 * t; + output[1, 1] = 2 * t * t; + output[2, 2] = 2; + output[2, 3] = 2 * Math.Sin(t); + output[3, 2] = 2 * Math.Sin(t); + output[3, 3] = 2 * Math.Pow(Math.Sin(t), 2); + } + + public override double ItemValue(Vector x, int itemIndex) + { + var ii = itemIndex + 1; + var t = ii / 5.0; + return Math.Pow(x[0] + t * x[1] - Math.Exp(t), 2.0) + Math.Pow(x[2] + x[3] * Math.Sin(t) - Math.Cos(t), 2); + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/BrownBadlyScaledFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/BrownBadlyScaledFunction.cs new file mode 100644 index 00000000..e0056244 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/BrownBadlyScaledFunction.cs @@ -0,0 +1,131 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class BrownBadlyScaledFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new BrownBadlyScaledFunction(), + InitialGuess = new double[] { 1, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 1e6, 2e-6 }, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new BrownBadlyScaledFunction(), + InitialGuess = new double[] { 1, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 1e6, 2e-6 }, + LowerBound = new double[] { -1e8, -1e8 }, + UpperBound = new double[] { 1e8, 1e8 }, + CaseName = "loose bounds" + }; + yield return new TestCase() + { + Function = new BrownBadlyScaledFunction(), + InitialGuess = new double[] { 1, 1 }, + MinimalValue = 0.784e3, + MinimizingPoint = new double[] { 1e6, 2e-6 }, + LowerBound = new double[] { 0, 3e-5 }, + UpperBound = new double[] { 1e6, 100 }, + CaseName = "tight bounds" + }; + } + } + + public BrownBadlyScaledFunction() { } + + public override string Description + { + get + { + return "Brown badly scaled fun (MGH #4)"; + } + } + + public override int ItemDimension + { + get + { + return 3; + } + } + + public override int ParameterDimension + { + get + { + return 2; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + switch (itemIndex) + { + case 0: + output[0] = 1; + output[1] = 0; + break; + case 1: + output[0] = 0; + output[1] = 1; + break; + case 2: + output[0] = x[1]; + output[1] = x[0]; + break; + default: + throw new ArgumentException("itemIndex must be <= 2"); + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + switch (itemIndex) + { + case 0: + case 1: + output[0, 0] = 0; + output[0, 1] = 0; + output[1, 0] = 0; + output[1, 1] = 0; + break; + case 2: + output[0, 0] = 0; + output[0, 1] = 1; + output[1, 0] = 1; + output[1, 1] = 0; + break; + default: + throw new ArgumentException("itemIndex must be <= 2"); + } + } + + public override double ItemValue(Vector x, int itemIndex) + { + switch (itemIndex) + { + case 0: + return x[0] - 1e6; + case 1: + return x[1] - 2e-6; + case 2: + return x[0] * x[1] - 2; + default: + throw new ArgumentException("itemIndex must be <= 2"); + } + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/FreudensteinAndRothFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/FreudensteinAndRothFunction.cs new file mode 100644 index 00000000..16cbe4c2 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/FreudensteinAndRothFunction.cs @@ -0,0 +1,98 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; +using DenseVector = MathNet.Numerics.LinearAlgebra.Double.DenseVector; +using DenseMatrix = MathNet.Numerics.LinearAlgebra.Double.DenseMatrix; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class FreudensteinAndRothFunction : BaseTestFunction + { + + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new FreudensteinAndRothFunction(), + InitialGuess = new double[] { 0.5, -2 }, + MinimizingPoint = new double[] { 5, 4 }, + MinimalValue = 0, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new FreudensteinAndRothFunction(), + InitialGuess = new double[] { 0.5, -2 }, + MinimizingPoint = new double[] {5, 4}, + MinimalValue = 0, + LowerBound = new double[] { -1000, -1000 }, + UpperBound = new double[] { 1000, 1000}, + CaseName = "loose bounds" + }; + } + } + + public override string Description { get { return "Freudenstein & Roth fun (MGH #2)"; } } + + public override int ParameterDimension + { + get + { + return 2; + } + } + + public override int ItemDimension + { + get + { + return 2; + } + } + + public override double ItemValue(Vector x, int itemIndex) + { + if (itemIndex == 0) + return -13 + x[0] + ((5 - x[1]) * x[1] - 2) * x[1]; + else + return -29 + x[0] + ((x[1] + 1) * x[1] - 14) * x[1]; + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + if (itemIndex == 0) + { + output[0] = 1; + output[1] = -2 + (5 - 2 * x[1]) * x[1] + (5 - x[1]) * x[1]; + } + else + { + output[0] = 1; + output[1] = -14 + x[1] * (1 + x[1]) + x[1] * (1 + 2 * x[1]); + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + if (itemIndex == 0) + { + output[0, 0] = 0; + output[0, 1] = 0; + output[1, 0] = 0; + output[1, 1] = 10 - 6 * x[1]; + } + else + { + output[0, 0] = 0; + output[0, 1] = 0; + output[1, 0] = 0; + output[1, 1] = 2 + 6 * x[1]; + } + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/HelicalValleyFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/HelicalValleyFunction.cs new file mode 100644 index 00000000..987eb7fe --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/HelicalValleyFunction.cs @@ -0,0 +1,178 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class HelicalValleyFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new HelicalValleyFunction(), + InitialGuess = new double[] { -1, 0, 0 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 1, 0, 0 }, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new HelicalValleyFunction(), + InitialGuess = new double[] { -1, 0, 0 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 1, 0, 0 }, + LowerBound = new double[] { -1000, -1000, -1000 }, + UpperBound = new double[] { 1000, 1000, 1000 }, + CaseName = "loose bounds" + }; + yield return new TestCase() + { + Function = new HelicalValleyFunction(), + InitialGuess = new double[] { -1, 0, 0 }, + MinimalValue = 0.99042212, + LowerBound = new double[] { -100, -1, -1 }, + UpperBound = new double[] { 0.8, 1, 1 }, + CaseName = "tight bounds" + }; + } + } + + public override string Description + { + get + { + return "Helical valley fun (MGH #7)"; + } + } + + public override int ItemDimension + { + get + { + return 3; + } + } + + public override int ParameterDimension + { + get + { + return 3; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + switch (itemIndex) + { + case 0: + output[0] = -100 * theta10(x[0], x[1]); + output[1] = -100 * theta01(x[0], x[1]); + output[2] = 10; + break; + case 1: + output[0] = (10 * x[0]) / Math.Sqrt(x[0]*x[0] + x[1]*x[1]); + output[1] = (10 * x[1]) / Math.Sqrt(x[0]*x[0] + x[1]*x[1]); + output[2] = 0; + break; + case 2: + output[0] = 0; + output[1] = 0; + output[2] = 1; + break; + default: + throw new ArgumentException("itemIndex must be <= 2"); + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + switch (itemIndex) + { + case 0: + output[0, 0] = -100 * theta20(x[0], x[1]); + output[0, 1] = -100 * theta11(x[0], x[1]); + output[0, 2] = 0; + output[1, 0] = -100 * theta11(x[0], x[1]); + output[1, 1] = -100 * theta02(x[0], x[1]); + output[1, 2] = 0; + output[2, 0] = 0; + output[2, 1] = 0; + output[2, 2] = 0; + break; + case 1: + output[0, 0] = (10 * x[1]*x[1]) / Math.Pow(x[0]*x[0] + x[1]*x[1],1.5); + output[0, 1] = (-10 * x[0] * x[1]) / Math.Pow(x[0]*x[0] + x[1]*x[1],1.5); + output[0, 2] = 0; + output[1, 0] = (-10 * x[0] * x[1]) / Math.Pow(x[0] * x[0] + x[1] * x[1], 1.5); + output[1, 1] = (10 * x[0]*x[0]) / Math.Pow(x[0] * x[0] + x[1] * x[1], 1.5); + output[1, 2] = 0; + output[2, 0] = 0; + output[2, 1] = 0; + output[2, 2] = 0; + break; + case 2: + for (int ii = 0; ii < 2; ++ii) + for (int jj = 0; jj < 2; ++jj) + output[ii, jj] = 0; + break; + default: + throw new ArgumentException("itemIndex must be <= 2"); + } + } + + private static double theta(double x1, double x2) + { + if (x1 >= 0) + return 0.5 * Math.Atan(x2 / x1) / Math.PI; + else + return 0.5 * Math.Atan(x2 / x1) / Math.PI + 0.5; + } + + private static double theta10(double x1, double x2) + { + return -(x2 / (2 * Math.PI * Math.Pow(x1,2) + 2 * Math.PI * Math.Pow(x2,2))); + } + + private static double theta01(double x1, double x2) + { + return x1 / (2 * Math.PI * x1*x1 + 2 * Math.PI * x2*x2); + } + + private static double theta20(double x1,double x2) + { + return (x1 * x2) / (Math.PI * Math.Pow(x1 * x1 + x2 * x2, 2)); + } + + private static double theta11(double x1, double x2) + { + return (-x1 * x1 + x2 * x2) / (2 * Math.PI * Math.Pow(x1 * x1 + x2 * x2, 2)); + } + + private static double theta02(double x1, double x2) + { + return -((x1 * x2) / (Math.PI * Math.Pow(x1*x1 + x2*x2, 2))); + } + + public override double ItemValue(Vector x, int itemIndex) + { + switch (itemIndex) + { + case 0: + return 10 * (x[2] - 10 * theta(x[0], x[1])); + case 1: + return 10 * (Math.Sqrt(x[0] * x[0] + x[1] * x[1]) - 1); + case 2: + return x[2]; + default: + throw new ArgumentException("itemIndex must be <= 2"); + } + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/ITestFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/ITestFunction.cs new file mode 100644 index 00000000..563faec2 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/ITestFunction.cs @@ -0,0 +1,69 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; +using DenseVector = MathNet.Numerics.LinearAlgebra.Double.DenseVector; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class TestCase + { + public string CaseName; + public ITestFunction Function; + public DenseVector InitialGuess; + public DenseVector LowerBound; + public DenseVector UpperBound; + public double MinimalValue; + public DenseVector MinimizingPoint; + + public bool IsBounded + { + get + { + return this.LowerBound != null && this.UpperBound != null; + } + } + + public bool IsUnbounded + { + get + { + return this.IsUnboundedOverride ?? this.LowerBound == null || this.UpperBound == null; + } + } + + public bool? IsUnboundedOverride; + + public string FullName + { + get + { + return $"{this.Function.Description} {this.CaseName}"; + } + } + } + + public interface ITestFunction + { + string Description { get; } + int ParameterDimension { get; } + int ItemDimension { get; } + + double ItemValue(Vector x, int itemIndex); + Vector ItemGradient(Vector x, int itemIndex); + void ItemGradientByRef(Vector x, int itemIndex, Vector output); + Matrix ItemHessian(Vector x, int itemIndex); + void ItemHessianByRef(Vector x, int itemIndex, Matrix output); + + Matrix Jacobian(Vector x); + void JacobianbyRef(Vector x, Matrix output); + + double SsqValue(Vector x); + Vector SsqGradient(Vector x); + void SsqGradientByRef(Vector x, Vector output); + Matrix SsqHessian(Vector x); + void SsqHessianByRef(Vector x, Matrix output); + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/JennrichAndSampsonFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/JennrichAndSampsonFunction.cs new file mode 100644 index 00000000..58cc30eb --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/JennrichAndSampsonFunction.cs @@ -0,0 +1,102 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class JennrichAndSampsonFunction : BaseTestFunction + { + private readonly int _m; + + public JennrichAndSampsonFunction(int itemDimension) + { + if (itemDimension < 2) + throw new ArgumentException("itemDimension must be at least 2."); + _m = itemDimension; + } + + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new JennrichAndSampsonFunction(10), + InitialGuess = new double[] { 0.3, 0.4 }, + MinimalValue = 124.362, + MinimizingPoint = new double[] { 0.2578, 0.2578 }, + CaseName = "unbounded" + }; + //yield return new TestCase() + //{ + // Function = new JennrichAndSampsonFunction(10), + // LowerBound = new double[] { 0.6, 0.5 }, + // UpperBound = new double[] { 10, 50 }, + // StartPoint = new double[] { 1.0, 1.0 }, + // MinimizingInput = null, + // MinimizingValue = 0, + // CaseName = "tight bounds" + //}; + yield return new TestCase() + { + Function = new JennrichAndSampsonFunction(10), + LowerBound = new double[] { -50, -50 }, + UpperBound = new double[] { 50, 50 }, + InitialGuess = new double[] { 0.3, 0.4 }, + MinimizingPoint = null, + MinimalValue = 0, + CaseName = "loose bounds" + }; + } + } + + public override string Description + { + get + { + return $"Jennrich & Sampson fun (MGH #6) (n={this.ItemDimension})"; + } + } + + public override int ItemDimension + { + get + { + return _m; + } + } + + public override int ParameterDimension + { + get + { + return 2; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + int ii = itemIndex + 1; + output[0] = -(Math.Exp(ii * x[0]) * ii); + output[1] = -(Math.Exp(ii * x[1]) * ii); + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + int ii = itemIndex + 1; + output[0, 0] = -(Math.Exp(ii * x[0]) * ii*ii); + output[0, 1] = 0; + output[1, 0] = 0; + output[1, 1] = -(Math.Exp(ii * x[1]) * ii*ii); + } + + public override double ItemValue(Vector x, int itemIndex) + { + int ii = itemIndex + 1; + return 2 + 2 * ii - (Math.Exp(ii * x[0]) + Math.Exp(ii * x[1])); + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/MeyerFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/MeyerFunction.cs new file mode 100644 index 00000000..961da226 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/MeyerFunction.cs @@ -0,0 +1,99 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class MeyerFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new MeyerFunction(), + InitialGuess = new double[] { 0.02, 4000, 250 }, + MinimalValue = 87.9458, + MinimizingPoint = null, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new MeyerFunction(), + InitialGuess = new double[] { 0.02, 4000, 250 }, + MinimalValue = 87.9458, + MinimizingPoint = null, + LowerBound = new double[] { -1e6, -1e6, -1e6 }, + UpperBound = new double[] { 1e6, 1e6, 1e6 }, + CaseName = "loose bounds" + }; + } + } + + public MeyerFunction() { } + + public override string Description + { + get + { + return "Meyer fun (MGH #10)"; + } + } + + public override int ItemDimension + { + get + { + return 16; + } + } + + public override int ParameterDimension + { + get + { + return 3; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + int ii = itemIndex + 1; + output[0] = Math.Exp(x[1] / (45.0 + 5 * ii + x[2])); + output[1] = (Math.Exp(x[1] / (45.0 + 5 * ii + x[2])) * x[0]) / (45 + 5 * ii + x[2]); + output[2] = -(Math.Exp(x[1] / (45.0 + 5 * ii + x[2])) * x[0] * x[1]) / Math.Pow(45 + 5 * ii + x[2], 2); + + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + var ii = itemIndex + 1; + + var t0 = (45.0 + 5 * ii + x[2]); + var t1 = Math.Exp(x[1] / t0); + + output[0, 0] = 0; + output[0, 1] = t1 / t0; + output[0, 2] = -t1 * x[1] / Math.Pow(t0, 2); + output[1, 0] = t1 / t0; + output[1, 1] = t1 * x[0] / Math.Pow(t0, 2); + output[1, 2] = -t1 * x[0] * (t0 + x[1]) / Math.Pow(t0, 3); + output[2, 0] = -t1 * x[1] / Math.Pow(t0, 2); + output[2, 1] = -t1 * x[0] * (t0 + x[1]) / Math.Pow(t0, 3); + output[2, 2] = t1 * x[0] * x[1] * (2*t0 + x[1]) / Math.Pow(t0, 4); + } + + private static readonly double[] y = { 34780, 28610, 23650, 19630, 16370, 13720, 11540, 9744, 8261, 7030, 6005, 5147, 4427, 3820, 3307, 2872 }; + + public override double ItemValue(Vector x, int itemIndex) + { + var ii = itemIndex + 1; + var t = 45.0 + 5 * ii; + return x[0] * Math.Exp(x[1] / (t + x[2])) - y[itemIndex]; + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/PowellBadlyScaledFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/PowellBadlyScaledFunction.cs new file mode 100644 index 00000000..8885a39d --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/PowellBadlyScaledFunction.cs @@ -0,0 +1,111 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class PowellBadlyScaledFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new PowellBadlyScaledFunction(), + InitialGuess = new double[] { 0, 1 }, + MinimizingPoint = new double[] { 1.098e-5, 9.106 }, + MinimalValue = 0, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new PowellBadlyScaledFunction(), + InitialGuess = new double[] { 0, 1 }, + MinimizingPoint = new double[] { 1.098e-5, 9.106 }, + MinimalValue = 0, + LowerBound = new double[] { -1000, -1000 }, + UpperBound = new double[] { 1000, 1000 }, + CaseName = "loose bounds" + }; + yield return new TestCase() + { + Function = new PowellBadlyScaledFunction(), + LowerBound = new double[] { 0, 1 }, + UpperBound = new double[] { 1, 9 }, + InitialGuess = new double[] { 0, 1 }, + MinimalValue = 0.15125900e-9, + CaseName = "tight bounds" + }; + } + } + + public override string Description + { + get + { + return "Powell badly scaled fun (MGH #3)"; + } + } + + public override int ItemDimension + { + get + { + return 2; + } + } + + public override int ParameterDimension + { + get + { + return 2; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + if (itemIndex == 0) + { + output[0] = 10000 * x[1]; + output[1] = 10000 * x[0]; + } + else if (itemIndex == 1) + { + output[0] = -Math.Exp(-x[0]); + output[1] = -Math.Exp(-x[1]); + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + if (itemIndex == 0) + { + output[0, 0] = 0; + output[0, 1] = 10000; + output[1, 0] = 10000; + output[1, 1] = 0; + } + else + { + output[0, 0] = Math.Exp(-x[0]); + output[0, 1] = 0; + output[1, 0] = 0; + output[1,1] = Math.Exp(-x[1]); + } + } + + public override double ItemValue(Vector x, int itemIndex) + { + if (itemIndex == 0) + return 10000.0 * x[0] * x[1] - 1; + else + return Math.Exp(-x[0]) + Math.Exp(-x[1]) - 1.0001; + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/PowellSingularFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/PowellSingularFunction.cs new file mode 100644 index 00000000..130b057a --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/PowellSingularFunction.cs @@ -0,0 +1,141 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class PowellSingularFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new PowellSingularFunction(), + InitialGuess = new double[] { 3, -1, 0, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] {0,0,0,0}, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new PowellSingularFunction(), + InitialGuess = new double[] { 3, -1, 0, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 0, 0, 0, 0 }, + LowerBound = new double[] {-1000, -1000, -1000, -1000}, + UpperBound = new double[] { 1000, 1000, 1000, 1000 }, + CaseName = "loose bounds" + }; + } + } + + + public PowellSingularFunction() { } + + public override string Description + { + get + { + return "Powell singular fun (MGH #13)"; + } + } + + public override int ItemDimension + { + get + { + return 4; + } + } + + public override int ParameterDimension + { + get + { + return 4; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + switch (itemIndex) + { + case 0: + output[0] = 1; + output[1] = 10; + output[2] = 0; + output[3] = 0; + break; + case 1: + output[0] = 0; + output[1] = 0; + output[2] = Math.Sqrt(5); + output[3] = -Math.Sqrt(5); + break; + case 2: + output[0] = 0; + output[1] = 2*(x[1]-2*x[2]); + output[2] = -4*x[1] + 8*x[2]; + output[3] = 0; + break; + case 3: + output[0] = 2*Math.Sqrt(10)*(x[0] - x[3]); + output[1] = 0; + output[2] = 0; + output[3] = -2*Math.Sqrt(10)*(x[0] - x[3]); + break; + default: + throw new ArgumentException("itemIndex must be <= 3"); + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + for (int ii = 0; ii < 4; ++ii) + for (int jj = 0; jj < 4; ++jj) + output[ii, jj] = 0; + switch(itemIndex) + { + case 0: + case 1: + break; + case 2: + output[1, 1] = 2; + output[1, 2] = -4; + output[2, 1] = -4; + output[2, 2] = 8; + break; + case 3: + output[0, 0] = 2 * Math.Sqrt(10); + output[0, 3] = -2 * Math.Sqrt(10); + output[3, 0] = -2 * Math.Sqrt(10); + output[3, 3] = 2 * Math.Sqrt(10); + break; + default: + throw new ArgumentException("itemIndex must be <= 3"); + } + } + + public override double ItemValue(Vector x, int itemIndex) + { + switch (itemIndex) + { + case 0: + return x[0] + 10 * x[1]; + case 1: + return Math.Sqrt(5) * (x[2] - x[3]); + case 2: + return Math.Pow(x[1] - 2 * x[2], 2); + case 3: + return Math.Sqrt(10.0) * Math.Pow(x[0] - x[3], 2); + default: + throw new ArgumentException("itemIndex must be <= 3"); + } + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/RosenbrockFunction2.cs b/src/UnitTests/OptimizationTests/TestFunctions/RosenbrockFunction2.cs new file mode 100644 index 00000000..c693b507 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/RosenbrockFunction2.cs @@ -0,0 +1,156 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class RosenbrockFunction2 : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { -1.2, 1 }, + MinimizingPoint = new double[] { 1, 1 }, + MinimalValue = 0, + LowerBound = new double[] { -1000, -1000 }, + UpperBound = new double[] { 1000, 1000 }, + CaseName = "hard start", + IsUnboundedOverride = true + }; + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { 1.2, 1.2 }, + MinimizingPoint = new double[] { 1, 1 }, + MinimalValue = 0, + LowerBound = new double[] { -5, -5 }, + UpperBound = new double[] { 5, 5 }, + CaseName = "easy start" + }; + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { -0.9, -0.5 }, + MinimizingPoint = new double[] { 1, 1 }, + MinimalValue = 0, + LowerBound = new double[] { -5, -5 }, + UpperBound = new double[] { 5, 5 }, + CaseName = "Overton start", + IsUnboundedOverride = true + }; + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { 1.2, 1.2 }, + MinimizingPoint = new double[] { 1, 1 }, + MinimalValue = 0, + LowerBound = new double[] { 1, -5 }, + UpperBound = new double[] { 5, 5 }, + CaseName = "easy one active bound" + }; + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { 1.2, 1.2 }, + MinimizingPoint = new double[] { 1, 1 }, + MinimalValue = 0, + LowerBound = new double[] { 1, 1 }, + UpperBound = new double[] { 5, 5 }, + CaseName = "easy two active bounds" + }; + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { 2.5, 2.5 }, + MinimizingPoint = new double[] { 2, 4 }, + MinimalValue = 1, + LowerBound = new double[] { 2, 2 }, + UpperBound = new double[] { 5, 5 }, + CaseName = "min on lower bound, not local" + }; + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { -0.9, -0.5 }, + MinimizingPoint = new double[] { 0.5, 0.25 }, + MinimalValue = 0.25, + LowerBound = new double[] { -2, -2 }, + UpperBound = new double[] { 0.5, 0.5 }, + CaseName = "min on upper bound, not local" + }; + + + } + } + public RosenbrockFunction2() { } + + public override string Description + { + get + { + return "Rosenbrock fun (MGH #1)"; + } + } + + public override int ItemDimension + { + get + { + return 2; + } + } + + public override int ParameterDimension + { + get + { + return 2; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + if (itemIndex == 0) + { + output[0] = -20 * x[0]; + output[1] = 10; + } else + { + output[0] = -1; + output[1] = 0; + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + if (itemIndex == 0) + { + output[0, 0] = -20; + output[0, 1] = 0; + output[1, 0] = 0; + output[1, 1] = 0; + } else + { + output[0, 0] = 0; + output[0, 1] = 0; + output[1, 0] = 0; + output[1, 1] = 0; + } + } + + public override double ItemValue(Vector x, int itemIndex) + { + if (itemIndex == 0) + return 10 * (x[1] - x[0] * x[0]); + else + return 1 - x[0]; + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/WoodFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/WoodFunction.cs new file mode 100644 index 00000000..4e8f014e --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/WoodFunction.cs @@ -0,0 +1,164 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class WoodFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new WoodFunction(), + InitialGuess = new double[] { -3, -1, -3, -1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 1, 1, 1, 1 }, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new WoodFunction(), + InitialGuess = new double[] { -3, -1, -3, -1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 1, 1, 1, 1 }, + LowerBound = new double[] { -1000, -1000, -1000, -1000 }, + UpperBound = new double[] { 1000, 1000, 1000, 1000 }, + CaseName = "loose bounds" + }; + yield return new TestCase() + { + Function = new WoodFunction(), + InitialGuess = new double[] { -3, -1, -3, -1 }, + MinimalValue = 1.5567008, + MinimizingPoint = null, + LowerBound = new double[] { -100, -100, -100, -100 }, + UpperBound = new double[] { 0, 10, 100, 100 }, + CaseName = "tight bounds" + }; + } + } + + public WoodFunction() { } + + public override string Description + { + get + { + return "Wood fun (MGH #14)"; + } + } + + public override int ItemDimension + { + get + { + return 6; + } + } + + public override int ParameterDimension + { + get + { + return 4; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + switch (itemIndex) + { + case 0: + output[0] = -20 * x[0]; + output[1] = 10; + output[2] = 0; + output[3] = 0; + break; + case 1: + output[0] = -1; + output[1] = 0; + output[2] = 0; + output[3] = 0; + break; + case 2: + output[0] = 0; + output[1] = 0; + output[2] = -6 * Math.Sqrt(10) * x[2]; + output[3] = 3 * Math.Sqrt(10); + break; + case 3: + output[0] = 0; + output[1] = 0; + output[2] = -1; + output[3] = 0; + break; + case 4: + output[0] = 0; + output[1] = Math.Sqrt(10); + output[2] = 0; + output[3] = Math.Sqrt(10); + break; + case 5: + output[0] = 0; + output[1] = 1.0 / Math.Sqrt(10); + output[2] = 0; + output[3] = -1.0 / Math.Sqrt(10); + break; + default: + throw new ArgumentException("itemIndex must be <= 5"); + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + for (int ii = 0; ii < 4; ++ii) + for (int jj = 0; jj < 4; ++jj) + output[ii, jj] = 0; + switch (itemIndex) + { + case 0: + output[0, 0] = -20; + break; + case 1: + break; + case 2: + output[2, 2] = -6 * Math.Sqrt(10); + break; + case 3: + case 4: + case 5: + break; + default: + throw new ArgumentException("itemIndex must be <= 5"); + + } + } + + public override double ItemValue(Vector x, int itemIndex) + { + switch (itemIndex) + { + case 0: + return 10 * (x[1] - x[0] * x[0]); + case 1: + return 1 - x[0]; + case 2: + return Math.Sqrt(90) * (x[3] - x[2] * x[2]); + case 3: + return 1 - x[2]; + case 4: + return Math.Sqrt(10) * (x[1] + x[3] - 2); + case 5: + return (x[1] - x[3]) / Math.Sqrt(10); + default: + throw new ArgumentException("itemIndex must be <= 5"); + } + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index cb59e924..8e3e3c98 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -345,16 +345,31 @@ + + - + + + + + + + + + + + + + + - - - - + + + + @@ -401,7 +416,7 @@ - +