From adfb2dc6b57eb6fa9b2f8246997c155a7206d9e4 Mon Sep 17 00:00:00 2001 From: Erik Ovegard Date: Sun, 30 Aug 2015 10:43:16 +0200 Subject: [PATCH 1/4] Initial conversion, tests remaining --- src/Numerics/Numerics.csproj | 1 + .../Optimization/MinimizationResult.cs | 3 +- .../Optimization/NelderMeadSimplex.cs | 339 ++++++++++++++++++ 3 files changed, 342 insertions(+), 1 deletion(-) create mode 100644 src/Numerics/Optimization/NelderMeadSimplex.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index eb3250c3..cae09415 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -102,6 +102,7 @@ + diff --git a/src/Numerics/Optimization/MinimizationResult.cs b/src/Numerics/Optimization/MinimizationResult.cs index 77cc389c..2a9d8fcd 100644 --- a/src/Numerics/Optimization/MinimizationResult.cs +++ b/src/Numerics/Optimization/MinimizationResult.cs @@ -13,7 +13,8 @@ namespace MathNet.Numerics.Optimization WeakWolfeCriteria, BoundTolerance, StrongWolfeCriteria, - LackOfFunctionImprovement + LackOfFunctionImprovement, + Converged } public Vector MinimizingPoint { get { return FunctionInfoAtMinimum.Point; } } diff --git a/src/Numerics/Optimization/NelderMeadSimplex.cs b/src/Numerics/Optimization/NelderMeadSimplex.cs new file mode 100644 index 00000000..583ce8de --- /dev/null +++ b/src/Numerics/Optimization/NelderMeadSimplex.cs @@ -0,0 +1,339 @@ +using MathNet.Numerics.LinearAlgebra; +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + public sealed class NelderMeadSimplex + { + private static readonly double JITTER = 1e-10d; // a small value used to protect against floating point noise + + public static MinimizationResult Regress(IObjectiveFunction objectiveFunction, Vector initialGuess, + double convergenceTolerance, int maxEvaluations) + { + SimplexConstant[] simplexConstants = SimplexConstant.CreateFromVector(initialGuess); + + // confirm that we are in a position to commence + if (objectiveFunction == null) + throw new InvalidOperationException("ObjectiveFunction must be set to a valid ObjectiveFunctionDelegate"); + + if (simplexConstants == null) + throw new InvalidOperationException("SimplexConstants must be initialized"); + + // create the initial simplex + int numDimensions = simplexConstants.Length; + int numVertices = numDimensions + 1; + Vector[] vertices = _initializeVertices(simplexConstants); + double[] errorValues = new double[numVertices]; + + int evaluationCount = 0; + MinimizationResult.ExitCondition exitCondition = MinimizationResult.ExitCondition.None; + ErrorProfile errorProfile; + + errorValues = _initializeErrorValues(vertices, objectiveFunction); + + // iterate until we converge, or complete our permitted number of iterations + while (true) + { + errorProfile = _evaluateSimplex(errorValues); + + // see if the range in point heights is small enough to exit + if (_hasConverged(convergenceTolerance, errorProfile, errorValues)) + { + exitCondition = MinimizationResult.ExitCondition.Converged; + break; + } + + // attempt a reflection of the simplex + double reflectionPointValue = _tryToScaleSimplex(-1.0, ref errorProfile, vertices, errorValues, objectiveFunction); + ++evaluationCount; + if (reflectionPointValue <= errorValues[errorProfile.LowestIndex]) + { + // it's better than the best point, so attempt an expansion of the simplex + double expansionPointValue = _tryToScaleSimplex(2.0, ref errorProfile, vertices, errorValues, objectiveFunction); + ++evaluationCount; + } + else if (reflectionPointValue >= errorValues[errorProfile.NextHighestIndex]) + { + // it would be worse than the second best point, so attempt a contraction to look + // for an intermediate point + double currentWorst = errorValues[errorProfile.HighestIndex]; + double contractionPointValue = _tryToScaleSimplex(0.5, ref errorProfile, vertices, errorValues, objectiveFunction); + ++evaluationCount; + if (contractionPointValue >= currentWorst) + { + // that would be even worse, so let's try to contract uniformly towards the low point; + // don't bother to update the error profile, we'll do it at the start of the + // next iteration + _shrinkSimplex(errorProfile, vertices, errorValues, objectiveFunction); + evaluationCount += numVertices; // that required one function evaluation for each vertex; keep track + } + } + // check to see if we have exceeded our alloted number of evaluations + if (evaluationCount >= maxEvaluations) + { + exitCondition = MinimizationResult.ExitCondition.LackOfProgress; + break; + } + } + var regressionResult = new MinimizationResult(null, evaluationCount, exitCondition); + return regressionResult; + } + + /// + /// Evaluate the objective function at each vertex to create a corresponding + /// list of error values for each vertex + /// + /// + /// + private static double[] _initializeErrorValues(Vector[] vertices, IObjectiveFunction objectiveFunction) + { + double[] errorValues = new double[vertices.Length]; + for (int i = 0; i < vertices.Length; i++) + { + objectiveFunction.EvaluateAt(vertices[i]); + errorValues[i] = objectiveFunction.Value; + } + return errorValues; + } + + /// + /// Check whether the points in the error profile have so little range that we + /// consider ourselves to have converged + /// + /// + /// + /// + private static bool _hasConverged(double convergenceTolerance, ErrorProfile errorProfile, double[] errorValues) + { + double range = 2 * Math.Abs(errorValues[errorProfile.HighestIndex] - errorValues[errorProfile.LowestIndex]) / + (Math.Abs(errorValues[errorProfile.HighestIndex]) + Math.Abs(errorValues[errorProfile.LowestIndex]) + JITTER); + + if (range < convergenceTolerance) + { + return true; + } + else + { + return false; + } + } + + /// + /// Examine all error values to determine the ErrorProfile + /// + /// + /// + private static ErrorProfile _evaluateSimplex(double[] errorValues) + { + ErrorProfile errorProfile = new ErrorProfile(); + if (errorValues[0] > errorValues[1]) + { + errorProfile.HighestIndex = 0; + errorProfile.NextHighestIndex = 1; + } + else + { + errorProfile.HighestIndex = 1; + errorProfile.NextHighestIndex = 0; + } + + for (int index = 0; index < errorValues.Length; index++) + { + double errorValue = errorValues[index]; + if (errorValue <= errorValues[errorProfile.LowestIndex]) + { + errorProfile.LowestIndex = index; + } + if (errorValue > errorValues[errorProfile.HighestIndex]) + { + errorProfile.NextHighestIndex = errorProfile.HighestIndex; // downgrade the current highest to next highest + errorProfile.HighestIndex = index; + } + else if (errorValue > errorValues[errorProfile.NextHighestIndex] && index != errorProfile.HighestIndex) + { + errorProfile.NextHighestIndex = index; + } + } + + return errorProfile; + } + + /// + /// Construct an initial simplex, given starting guesses for the constants, and + /// initial step sizes for each dimension + /// + /// + /// + private static Vector[] _initializeVertices(SimplexConstant[] simplexConstants) + { + int numDimensions = simplexConstants.Length; + Vector[] vertices = new Vector[numDimensions + 1]; + + // define one point of the simplex as the given initial guesses + var p0 = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(numDimensions); + for (int i = 0; i < numDimensions; i++) + { + p0[i] = simplexConstants[i].Value; + } + + // now fill in the vertices, creating the additional points as: + // P(i) = P(0) + Scale(i) * UnitVector(i) + vertices[0] = p0; + for (int i = 0; i < numDimensions; i++) + { + double scale = simplexConstants[i].InitialPerturbation; + Vector unitVector = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(numDimensions); + unitVector[i] = 1; + vertices[i + 1] = p0.Add(unitVector.Multiply(scale)); + } + return vertices; + } + + /// + /// Test a scaling operation of the high point, and replace it if it is an improvement + /// + /// + /// + /// + /// + /// + private static double _tryToScaleSimplex(double scaleFactor, ref ErrorProfile errorProfile, Vector[] vertices, + double[] errorValues, IObjectiveFunction objectiveFunction) + { + // find the centroid through which we will reflect + Vector centroid = _computeCentroid(vertices, errorProfile); + + // define the vector from the centroid to the high point + Vector centroidToHighPoint = vertices[errorProfile.HighestIndex].Subtract(centroid); + + // scale and position the vector to determine the new trial point + Vector newPoint = centroidToHighPoint.Multiply(scaleFactor).Add(centroid); + + // evaluate the new point + objectiveFunction.EvaluateAt(newPoint); + double newErrorValue = objectiveFunction.Value; + + // if it's better, replace the old high point + if (newErrorValue < errorValues[errorProfile.HighestIndex]) + { + vertices[errorProfile.HighestIndex] = newPoint; + errorValues[errorProfile.HighestIndex] = newErrorValue; + } + + return newErrorValue; + } + + /// + /// Contract the simplex uniformly around the lowest point + /// + /// + /// + /// + private static void _shrinkSimplex(ErrorProfile errorProfile, Vector[] vertices, double[] errorValues, + IObjectiveFunction objectiveFunction) + { + Vector lowestVertex = vertices[errorProfile.LowestIndex]; + for (int i = 0; i < vertices.Length; i++) + { + if (i != errorProfile.LowestIndex) + { + vertices[i] = (vertices[i].Add(lowestVertex)).Multiply(0.5); + objectiveFunction.EvaluateAt(vertices[i]); + errorValues[i] = objectiveFunction.Value; + } + } + } + + /// + /// Compute the centroid of all points except the worst + /// + /// + /// + /// + private static Vector _computeCentroid(Vector[] vertices, ErrorProfile errorProfile) + { + int numVertices = vertices.Length; + // find the centroid of all points except the worst one + Vector centroid = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(numVertices - 1); + for (int i = 0; i < numVertices; i++) + { + if (i != errorProfile.HighestIndex) + { + centroid = centroid.Add(vertices[i]); + } + } + return centroid.Multiply(1.0d / (numVertices - 1)); + } + + private sealed class SimplexConstant + { + private double _value; + private double _initialPerturbation; + + public SimplexConstant(double value, double initialPerturbation) + { + _value = value; + _initialPerturbation = initialPerturbation; + } + + /// + /// The value of the constant + /// + public double Value + { + get { return _value; } + set { _value = value; } + } + + // The size of the initial perturbation + public double InitialPerturbation + { + get { return _initialPerturbation; } + set { _initialPerturbation = value; } + } + + public static SimplexConstant[] CreateFromVector(Vector initialGuess) + { + var constants = new SimplexConstant[initialGuess.Count]; + + for (int i = 0; i < constants.Length;i++ ) + { + double pertubation = initialGuess[i]==0.0 ? 1e-5 : initialGuess[i]*1e-5; + constants[i] = new SimplexConstant(initialGuess[i], pertubation); + } + return constants; + } + } + + private sealed class ErrorProfile + { + private int _highestIndex; + private int _nextHighestIndex; + private int _lowestIndex; + + public int HighestIndex + { + get { return _highestIndex; } + set { _highestIndex = value; } + } + + public int NextHighestIndex + { + get { return _nextHighestIndex; } + set { _nextHighestIndex = value; } + } + + public int LowestIndex + { + get { return _lowestIndex; } + set { _lowestIndex = value; } + } + } + } + + + +} From 10e2a898c622903fbe7114a2cd97808d88cbe7c5 Mon Sep 17 00:00:00 2001 From: Erik Ovegard Date: Sun, 30 Aug 2015 13:07:54 +0200 Subject: [PATCH 2/4] Added tests and copyright notice --- .../Optimization/NelderMeadSimplex.cs | 66 +++++++++-- .../NelderMeadSimplexTests.cs | 107 ++++++++++++++++++ src/UnitTests/UnitTests.csproj | 1 + 3 files changed, 164 insertions(+), 10 deletions(-) create mode 100644 src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs diff --git a/src/Numerics/Optimization/NelderMeadSimplex.cs b/src/Numerics/Optimization/NelderMeadSimplex.cs index 583ce8de..2afa1113 100644 --- a/src/Numerics/Optimization/NelderMeadSimplex.cs +++ b/src/Numerics/Optimization/NelderMeadSimplex.cs @@ -1,4 +1,36 @@ -using MathNet.Numerics.LinearAlgebra; +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2015 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +// Converted from code relased with a MIT liscense available at https://code.google.com/p/nelder-mead-simplex/ + +using MathNet.Numerics.LinearAlgebra; using System; using System.Collections.Generic; using System.Linq; @@ -10,17 +42,31 @@ namespace MathNet.Numerics.Optimization { private static readonly double JITTER = 1e-10d; // a small value used to protect against floating point noise - public static MinimizationResult Regress(IObjectiveFunction objectiveFunction, Vector initialGuess, - double convergenceTolerance, int maxEvaluations) + public double ConvergenceTolerance { get; set; } + public int MaximumIterations { get; set; } + + public NelderMeadSimplex(double convergenceTolerance, int maximumIterations) { - SimplexConstant[] simplexConstants = SimplexConstant.CreateFromVector(initialGuess); + ConvergenceTolerance = convergenceTolerance; + MaximumIterations = maximumIterations; + } + /// + /// Finds the minimum of the objective function + /// + /// The objective function, no gradient or hessian needed + /// The intial guess + /// The minimum point + public MinimizationResult FindMinimum(IObjectiveFunction objectiveFunction, Vector initialGuess) + { // confirm that we are in a position to commence if (objectiveFunction == null) - throw new InvalidOperationException("ObjectiveFunction must be set to a valid ObjectiveFunctionDelegate"); + throw new ArgumentNullException("objectiveFunction","ObjectiveFunction must be set to a valid ObjectiveFunctionDelegate"); - if (simplexConstants == null) - throw new InvalidOperationException("SimplexConstants must be initialized"); + if (initialGuess == null) + throw new ArgumentNullException("initialGuess", "initialGuess must be initialized"); + + SimplexConstant[] simplexConstants = SimplexConstant.CreateFromVector(initialGuess); // create the initial simplex int numDimensions = simplexConstants.Length; @@ -40,7 +86,7 @@ namespace MathNet.Numerics.Optimization errorProfile = _evaluateSimplex(errorValues); // see if the range in point heights is small enough to exit - if (_hasConverged(convergenceTolerance, errorProfile, errorValues)) + if (_hasConverged(ConvergenceTolerance, errorProfile, errorValues)) { exitCondition = MinimizationResult.ExitCondition.Converged; break; @@ -72,13 +118,13 @@ namespace MathNet.Numerics.Optimization } } // check to see if we have exceeded our alloted number of evaluations - if (evaluationCount >= maxEvaluations) + if (evaluationCount >= MaximumIterations) { exitCondition = MinimizationResult.ExitCondition.LackOfProgress; break; } } - var regressionResult = new MinimizationResult(null, evaluationCount, exitCondition); + var regressionResult = new MinimizationResult(objectiveFunction, evaluationCount, exitCondition); return regressionResult; } diff --git a/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs b/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs new file mode 100644 index 00000000..f1a2e493 --- /dev/null +++ b/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs @@ -0,0 +1,107 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2015 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using NUnit.Framework; +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.Optimization; +using MathNet.Numerics.LinearAlgebra.Double; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + public class NelderMeadSimplexTests + { + /// + /// Test that finds the constants of a parable, function adds noise and return the mean square error + /// Copied from the test in https://code.google.com/p/nelder-mead-simplex/ + /// + [Test] + public void FindParableConstantsThatMinimizesErrors() + { + var nms = new NelderMeadSimplex(1e-6, 1000); + double a = 5; + double b = 10; + IObjectiveFunction objFun = ObjectiveFunction.Value((constants)=> + { + double ssq = 0; + System.Random r = new System.Random(); + for (double x = -10; x < 10; x += .1) + { + double yTrue = a * x * x + b * x + r.NextDouble(); + double yRegress = constants[0] * x * x + constants[1] * x; + ssq += Math.Pow((yTrue - yRegress), 2); + } + return ssq; + }); + var initialGuess = new DenseVector(2); + initialGuess[0] = 3; + initialGuess[1] = 5; + var result = nms.FindMinimum(objFun, initialGuess); + + Assert.NotNull(result); + Assert.NotNull(result.MinimizingPoint); + Assert.NotNull(result.FunctionInfoAtMinimum); + Assert.That(Math.Abs(result.MinimizingPoint[0] - a), Is.LessThan(1e-2)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - b), Is.LessThan(1e-2)); + } + + [Test] + public void NMS_FindMinimum_Rosenbrock_Easy() + { + var obj = ObjectiveFunction.Value(RosenbrockFunction.Value); + var solver = new NelderMeadSimplex(1e-5, maximumIterations: 1000); + var initialGuess = new DenseVector(new[] { 1.2, 1.2 }); + + var result = solver.FindMinimum(obj, initialGuess); + + 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 NMS_FindMinimum_Rosenbrock_Hard() + { + var obj = ObjectiveFunction.Value(RosenbrockFunction.Value); + var solver = new NelderMeadSimplex(1e-5, maximumIterations: 1000); + + var initialGuess = new DenseVector(new[] { -1.2, 1.0 }); + + var result = solver.FindMinimum(obj,initialGuess); + + 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/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 23b6164f..cb59e924 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -345,6 +345,7 @@ + From c942f60eef267d92794083454e3cfb27a46ada4d Mon Sep 17 00:00:00 2001 From: Erik Ovegard Date: Sun, 30 Aug 2015 13:35:57 +0200 Subject: [PATCH 3/4] Added more documentation to NelderMeadSimplex and changed style to match existing code --- .../Optimization/NelderMeadSimplex.cs | 74 +++++++++++++------ 1 file changed, 50 insertions(+), 24 deletions(-) diff --git a/src/Numerics/Optimization/NelderMeadSimplex.cs b/src/Numerics/Optimization/NelderMeadSimplex.cs index 2afa1113..72bdb2ed 100644 --- a/src/Numerics/Optimization/NelderMeadSimplex.cs +++ b/src/Numerics/Optimization/NelderMeadSimplex.cs @@ -38,6 +38,13 @@ using System.Text; namespace MathNet.Numerics.Optimization { + /// + /// Class implementing the Nelder-Mead simplex algorithm, used to find a minima when no gradient is available. + /// Called fminsearch() in Matlab. A description of the algorithm can be found at + /// http://se.mathworks.com/help/matlab/math/optimizing-nonlinear-functions.html#bsgpq6p-11 + /// or + /// https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method + /// public sealed class NelderMeadSimplex { private static readonly double JITTER = 1e-10d; // a small value used to protect against floating point noise @@ -52,12 +59,31 @@ namespace MathNet.Numerics.Optimization } /// - /// Finds the minimum of the objective function + /// Finds the minimum of the objective function without an intial pertubation, the default values used + /// by fminsearch() in Matlab are used instead + /// http://se.mathworks.com/help/matlab/math/optimizing-nonlinear-functions.html#bsgpq6p-11 /// /// The objective function, no gradient or hessian needed /// The intial guess /// The minimum point public MinimizationResult FindMinimum(IObjectiveFunction objectiveFunction, Vector initialGuess) + { + var initalPertubation = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(initialGuess.Count); + for (int i = 0; i < initialGuess.Count; i++) + { + initalPertubation[i] = initialGuess[i] == 0.0 ? 0.00025 : initialGuess[i] * 0.05; + } + return FindMinimum(objectiveFunction, initialGuess, initalPertubation); + } + + /// + /// Finds the minimum of the objective function with an intial pertubation + /// + /// The objective function, no gradient or hessian needed + /// The intial guess + /// The inital pertubation + /// The minimum point + public MinimizationResult FindMinimum(IObjectiveFunction objectiveFunction, Vector initialGuess, Vector initalPertubation) { // confirm that we are in a position to commence if (objectiveFunction == null) @@ -66,39 +92,42 @@ namespace MathNet.Numerics.Optimization if (initialGuess == null) throw new ArgumentNullException("initialGuess", "initialGuess must be initialized"); - SimplexConstant[] simplexConstants = SimplexConstant.CreateFromVector(initialGuess); + if (initialGuess == null) + throw new ArgumentNullException("initalPertubation", "initalPertubation must be initialized, if unknown use overloaded version of FindMinimum()"); + + SimplexConstant[] simplexConstants = SimplexConstant.CreateSimplexConstantsFromVectors(initialGuess,initalPertubation); // create the initial simplex int numDimensions = simplexConstants.Length; int numVertices = numDimensions + 1; - Vector[] vertices = _initializeVertices(simplexConstants); + Vector[] vertices = InitializeVertices(simplexConstants); double[] errorValues = new double[numVertices]; int evaluationCount = 0; MinimizationResult.ExitCondition exitCondition = MinimizationResult.ExitCondition.None; ErrorProfile errorProfile; - errorValues = _initializeErrorValues(vertices, objectiveFunction); + errorValues = InitializeErrorValues(vertices, objectiveFunction); // iterate until we converge, or complete our permitted number of iterations while (true) { - errorProfile = _evaluateSimplex(errorValues); + errorProfile = EvaluateSimplex(errorValues); // see if the range in point heights is small enough to exit - if (_hasConverged(ConvergenceTolerance, errorProfile, errorValues)) + if (HasConverged(ConvergenceTolerance, errorProfile, errorValues)) { exitCondition = MinimizationResult.ExitCondition.Converged; break; } // attempt a reflection of the simplex - double reflectionPointValue = _tryToScaleSimplex(-1.0, ref errorProfile, vertices, errorValues, objectiveFunction); + double reflectionPointValue = TryToScaleSimplex(-1.0, ref errorProfile, vertices, errorValues, objectiveFunction); ++evaluationCount; if (reflectionPointValue <= errorValues[errorProfile.LowestIndex]) { // it's better than the best point, so attempt an expansion of the simplex - double expansionPointValue = _tryToScaleSimplex(2.0, ref errorProfile, vertices, errorValues, objectiveFunction); + double expansionPointValue = TryToScaleSimplex(2.0, ref errorProfile, vertices, errorValues, objectiveFunction); ++evaluationCount; } else if (reflectionPointValue >= errorValues[errorProfile.NextHighestIndex]) @@ -106,22 +135,21 @@ namespace MathNet.Numerics.Optimization // it would be worse than the second best point, so attempt a contraction to look // for an intermediate point double currentWorst = errorValues[errorProfile.HighestIndex]; - double contractionPointValue = _tryToScaleSimplex(0.5, ref errorProfile, vertices, errorValues, objectiveFunction); + double contractionPointValue = TryToScaleSimplex(0.5, ref errorProfile, vertices, errorValues, objectiveFunction); ++evaluationCount; if (contractionPointValue >= currentWorst) { // that would be even worse, so let's try to contract uniformly towards the low point; // don't bother to update the error profile, we'll do it at the start of the // next iteration - _shrinkSimplex(errorProfile, vertices, errorValues, objectiveFunction); + ShrinkSimplex(errorProfile, vertices, errorValues, objectiveFunction); evaluationCount += numVertices; // that required one function evaluation for each vertex; keep track } } // check to see if we have exceeded our alloted number of evaluations if (evaluationCount >= MaximumIterations) { - exitCondition = MinimizationResult.ExitCondition.LackOfProgress; - break; + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); } } var regressionResult = new MinimizationResult(objectiveFunction, evaluationCount, exitCondition); @@ -134,7 +162,7 @@ namespace MathNet.Numerics.Optimization /// /// /// - private static double[] _initializeErrorValues(Vector[] vertices, IObjectiveFunction objectiveFunction) + private static double[] InitializeErrorValues(Vector[] vertices, IObjectiveFunction objectiveFunction) { double[] errorValues = new double[vertices.Length]; for (int i = 0; i < vertices.Length; i++) @@ -152,7 +180,7 @@ namespace MathNet.Numerics.Optimization /// /// /// - private static bool _hasConverged(double convergenceTolerance, ErrorProfile errorProfile, double[] errorValues) + private static bool HasConverged(double convergenceTolerance, ErrorProfile errorProfile, double[] errorValues) { double range = 2 * Math.Abs(errorValues[errorProfile.HighestIndex] - errorValues[errorProfile.LowestIndex]) / (Math.Abs(errorValues[errorProfile.HighestIndex]) + Math.Abs(errorValues[errorProfile.LowestIndex]) + JITTER); @@ -172,7 +200,7 @@ namespace MathNet.Numerics.Optimization /// /// /// - private static ErrorProfile _evaluateSimplex(double[] errorValues) + private static ErrorProfile EvaluateSimplex(double[] errorValues) { ErrorProfile errorProfile = new ErrorProfile(); if (errorValues[0] > errorValues[1]) @@ -213,7 +241,7 @@ namespace MathNet.Numerics.Optimization /// /// /// - private static Vector[] _initializeVertices(SimplexConstant[] simplexConstants) + private static Vector[] InitializeVertices(SimplexConstant[] simplexConstants) { int numDimensions = simplexConstants.Length; Vector[] vertices = new Vector[numDimensions + 1]; @@ -246,11 +274,11 @@ namespace MathNet.Numerics.Optimization /// /// /// - private static double _tryToScaleSimplex(double scaleFactor, ref ErrorProfile errorProfile, Vector[] vertices, + private static double TryToScaleSimplex(double scaleFactor, ref ErrorProfile errorProfile, Vector[] vertices, double[] errorValues, IObjectiveFunction objectiveFunction) { // find the centroid through which we will reflect - Vector centroid = _computeCentroid(vertices, errorProfile); + Vector centroid = ComputeCentroid(vertices, errorProfile); // define the vector from the centroid to the high point Vector centroidToHighPoint = vertices[errorProfile.HighestIndex].Subtract(centroid); @@ -278,7 +306,7 @@ namespace MathNet.Numerics.Optimization /// /// /// - private static void _shrinkSimplex(ErrorProfile errorProfile, Vector[] vertices, double[] errorValues, + private static void ShrinkSimplex(ErrorProfile errorProfile, Vector[] vertices, double[] errorValues, IObjectiveFunction objectiveFunction) { Vector lowestVertex = vertices[errorProfile.LowestIndex]; @@ -299,7 +327,7 @@ namespace MathNet.Numerics.Optimization /// /// /// - private static Vector _computeCentroid(Vector[] vertices, ErrorProfile errorProfile) + private static Vector ComputeCentroid(Vector[] vertices, ErrorProfile errorProfile) { int numVertices = vertices.Length; // find the centroid of all points except the worst one @@ -341,14 +369,12 @@ namespace MathNet.Numerics.Optimization set { _initialPerturbation = value; } } - public static SimplexConstant[] CreateFromVector(Vector initialGuess) + public static SimplexConstant[] CreateSimplexConstantsFromVectors(Vector initialGuess, Vector initialPertubation) { var constants = new SimplexConstant[initialGuess.Count]; - for (int i = 0; i < constants.Length;i++ ) { - double pertubation = initialGuess[i]==0.0 ? 1e-5 : initialGuess[i]*1e-5; - constants[i] = new SimplexConstant(initialGuess[i], pertubation); + constants[i] = new SimplexConstant(initialGuess[i], initialPertubation[i]); } return constants; } From ac51181682f5df917593d927f560315f0042d02c Mon Sep 17 00:00:00 2001 From: Erik Ovegard Date: Sun, 30 Aug 2015 14:03:46 +0200 Subject: [PATCH 4/4] Fixed XML-errors found by AppVeyor --- src/Numerics/Optimization/NelderMeadSimplex.cs | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/Numerics/Optimization/NelderMeadSimplex.cs b/src/Numerics/Optimization/NelderMeadSimplex.cs index 72bdb2ed..0024e96f 100644 --- a/src/Numerics/Optimization/NelderMeadSimplex.cs +++ b/src/Numerics/Optimization/NelderMeadSimplex.cs @@ -161,6 +161,7 @@ namespace MathNet.Numerics.Optimization /// list of error values for each vertex /// /// + /// /// private static double[] InitializeErrorValues(Vector[] vertices, IObjectiveFunction objectiveFunction) { @@ -177,6 +178,7 @@ namespace MathNet.Numerics.Optimization /// Check whether the points in the error profile have so little range that we /// consider ourselves to have converged /// + /// /// /// /// @@ -273,6 +275,7 @@ namespace MathNet.Numerics.Optimization /// /// /// + /// /// private static double TryToScaleSimplex(double scaleFactor, ref ErrorProfile errorProfile, Vector[] vertices, double[] errorValues, IObjectiveFunction objectiveFunction) @@ -306,6 +309,7 @@ namespace MathNet.Numerics.Optimization /// /// /// + /// private static void ShrinkSimplex(ErrorProfile errorProfile, Vector[] vertices, double[] errorValues, IObjectiveFunction objectiveFunction) {