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 @@ - +