diff --git a/src/Numerics/Optimization/BrentMinimizer.cs b/src/Numerics/Optimization/BrentMinimizer.cs
index 6b6de901..3d029903 100644
--- a/src/Numerics/Optimization/BrentMinimizer.cs
+++ b/src/Numerics/Optimization/BrentMinimizer.cs
@@ -5,7 +5,301 @@ using System.Text;
namespace MathNet.Numerics.Optimization
{
+ ///
+ /// Options for Brent Minimization.
+ ///
+ public class BrentOptions
+ {
+ public int MaximumIterations = 1000;
+
+ public double FunctionTolerance = 1e-4;
+ }
+
+ ///
+ /// Minimizes f(p) where p is a model parameter scalar, i.e. a line-search.
+ ///
public class BrentMinimizer
{
+ const double verySmallNumber = 1e-21, goldenRatio = 1.618034, minimumTolerance = 1.0e-11;
+ const double growLimit = 110.0, conjugateGradient = 0.3819660;
+ const int maxIterations = 500;
+
+ int bracketFunctionCalls;
+ public int FunctionCalls, Iterations;
+ double minPoint, minFunction;
+ Func function;
+ double pointA, pointB, pointC;
+ double functionA, functionB, functionC;
+
+ public double Tolerance { get; set; }
+
+ public BrentMinimizer(Func function)
+ {
+ this.function = function;
+ Tolerance = 1e-4;
+ }
+
+ public int Search(out double minPoint, out double minFunction)
+ {
+ bracketFunctionCalls = 0;
+ UpdateBracketInterval(0, 1);
+ FunctionCalls = 0;
+ Iterations = 0;
+ int result = BrentMinimize();
+ FunctionCalls += bracketFunctionCalls;
+ minPoint = this.minPoint;
+ minFunction = this.minFunction;
+ return result;
+ }
+
+ private int UpdateBracketInterval(double pointAStart, double pointBStart)
+ {
+ int iterations = 0;
+ pointA = pointAStart;
+ pointB = pointBStart;
+ int maxIterations;
+ maxIterations = 1000;
+ functionA = function(pointA);
+ functionB = function(pointB);
+ double temp;
+ if (functionA < functionB) // Swap points over
+ {
+ temp = functionA; functionA = functionB; functionB = temp;
+ temp = pointA; pointA = pointB; pointB = temp;
+ }
+ pointC = pointB + goldenRatio * (pointB - pointA);
+ functionC = function(pointC);
+ bracketFunctionCalls = 3; iterations = 0;
+ double temp1, temp2, value, denom;
+ while (functionC < functionB)
+ {
+ double pointW, functionW, wlim;
+ temp1 = (pointB - pointA) * (functionB - functionC);
+ temp2 = (pointB - pointC) * (functionB - functionA);
+ value = temp2 - temp1;
+ if (Math.Abs(value) < verySmallNumber) denom = 2.0 * verySmallNumber;
+ else denom = 2.0 * value;
+ pointW = pointB - ((pointB - pointC) * temp2 - (pointB - pointA) * temp1) / denom;
+ wlim = pointB + growLimit * (pointC - pointB);
+ if (iterations > maxIterations) return 1;
+ iterations++;
+ if ((pointW - pointC) * (pointB - pointW) > 0.0)
+ {
+ functionW = function(pointW);
+ bracketFunctionCalls++;
+ if (functionW < functionC)
+ {
+ pointA = pointB; pointB = pointW;
+ functionA = functionB; functionB = functionW;
+ return 0;
+ }
+ else if (functionW > functionB)
+ {
+ pointC = pointW; functionC = functionW;
+ return 0;
+ }
+ pointW = pointC + goldenRatio * (pointC - pointB);
+ functionW = function(pointW);
+ bracketFunctionCalls++;
+ }
+ else if ((pointW - wlim) * (wlim - pointC) >= 0.0)
+ {
+ pointW = wlim;
+ functionW = function(pointW);
+ bracketFunctionCalls++;
+ }
+ else if ((pointW - wlim) * (pointC - pointW) > 0.0)
+ {
+ functionW = function(pointW);
+ bracketFunctionCalls++;
+ if (functionW < functionC)
+ {
+ pointB = pointC; pointC = pointW;
+ pointW = pointC + goldenRatio * (pointC - pointB);
+ functionB = functionC; functionC = functionW;
+ functionW = function(pointW);
+ bracketFunctionCalls++;
+ }
+ }
+ else
+ {
+ pointW = pointC + goldenRatio * (pointC - pointB);
+ functionW = function(pointW);
+ bracketFunctionCalls++;
+ }
+ pointA = pointB; pointB = pointC; pointC = pointW;
+ functionA = functionB; functionB = functionC; functionC = functionW;
+ }
+ return 0;
+ }
+
+ // Find the minimum of the function using the Brent method with the current
+ // bracketing interval.
+ private int BrentMinimize()
+ {
+ int result = 0;
+ double x, w, v, fx, fw, fv;
+ double a, b, deltax, rat;
+ double cg = conjugateGradient;
+ x = w = v = pointB; // x is the point with lowest function value encountered
+ fw = fv = fx = function(x);
+ if (pointA < pointC)
+ {
+ a = pointA; b = pointC;
+ }
+ else
+ {
+ a = pointC; b = pointA;
+ }
+ deltax = 0.0;
+ FunctionCalls = 1;
+ Iterations = 0;
+ rat = 0;
+ while (Iterations < maxIterations)
+ {
+ double tol1, tol2, xmin, fval, xmid;
+ double temp1, temp2, p;
+ double u, fu, dx_temp;
+ tol1 = Tolerance * Math.Abs(x) + minimumTolerance;
+ tol2 = 2.0 * tol1;
+ xmid = 0.5 * (a + b);
+ if (Math.Abs(x - xmid) < (tol2 - 0.5 * (b - a))) // check for convergence
+ {
+ xmin = x; fval = fx;
+ result = 1;
+ break;
+ }
+ if (Math.Abs(deltax) <= tol1)
+ {
+ if (x >= xmid) deltax = a - x; // do a golden section step
+ else deltax = b - x;
+ rat = cg * deltax;
+ }
+ else // do a parabolic step
+ {
+ temp1 = (x - w) * (fx - fv);
+ temp2 = (x - v) * (fx - fw);
+ p = (x - v) * temp2 - (x - w) * temp1;
+ temp2 = 2.0 * (temp2 - temp1);
+ if (temp2 > 0.0) p = -p;
+ temp2 = Math.Abs(temp2);
+ dx_temp = deltax;
+ deltax = rat;
+ // check parabolic fit
+ if ((p > temp2 * (a - x)) && (p < temp2 * (b - x))
+ && (Math.Abs(p) < Math.Abs(0.5 * temp2 * dx_temp)))
+ {
+ rat = p * 1.0 / temp2; // if parabolic step is useful.
+ u = x + rat;
+ if (((u - a) < tol2) || ((b - u) < tol2))
+ {
+ if ((xmid - x) >= 0) rat = tol1;
+ else rat = -tol1;
+ }
+ }
+ else
+ {
+ if (x >= xmid) deltax = a - x; // if it's not do a golden section step
+ else deltax = b - x;
+ rat = cg * deltax;
+ }
+ }
+
+ if (Math.Abs(rat) < tol1) // update by at least tol1
+ {
+ if (rat >= 0) u = x + tol1;
+ else u = x - tol1;
+ }
+ else
+ {
+ u = x + rat;
+ }
+ fu = function(u);
+ FunctionCalls++;
+ if (fu > fx) // if it's bigger than current
+ {
+ if (u < x) a = u;
+ else b = u;
+ if ((fu <= fw) || (w == x))
+ {
+ v = w; w = u; fv = fw; fw = fu;
+ }
+ else if ((fu <= fv) || (v == x) || (v == w))
+ {
+ v = u; fv = fu;
+ }
+ }
+ else
+ {
+ if (u >= x) a = x;
+ else b = x;
+ v = w; w = x; x = u;
+ fv = fw; fw = fx; fx = fu;
+ }
+ Iterations++;
+ }
+ this.minPoint = x;
+ this.minFunction = fx;
+ return result;
+ }
+ }
+
+ ///
+ /// Minimizes f(p u) where p is a model parameter scalar and u is a direction vector.
+ ///
+ public class MultiDimensionalBrent
+ {
+ private Func function;
+ private int functionCalls;
+ double[] point;
+ BrentMinimizer lineSearch;
+
+ public MultiDimensionalBrent(Func function)
+ {
+ this.function = function;
+ lineSearch = new BrentMinimizer(this.PointAlongLine);
+ functionCalls = 0;
+ }
+
+ public double Tolerance
+ {
+ get { return lineSearch.Tolerance; }
+ set { lineSearch.Tolerance = value; }
+ }
+
+ public int FunctionCalls
+ {
+ get { return lineSearch.FunctionCalls; }
+ }
+
+ public double[] StartingPoint { get; set; }
+ public double[] Direction { get; set; }
+
+
+ public void SetDimension(int N)
+ {
+ point = new double[N];
+ }
+
+ public int Search(out double[] minPoint, out double minFunction)
+ {
+ double path;
+ int result = lineSearch.Search(out path, out minFunction);
+ for (int i = 0; i < StartingPoint.Length; ++i)
+ {
+ point[i] = StartingPoint[i] + path * Direction[i];
+ }
+ minPoint = point;
+ return result;
+ }
+
+ public double PointAlongLine(double path)
+ {
+ for (int i = 0; i < StartingPoint.Length; ++i)
+ {
+ point[i] = StartingPoint[i] + path * Direction[i];
+ }
+ return function(point);
+ }
}
}
diff --git a/src/Numerics/Optimization/PowellMinimizer.cs b/src/Numerics/Optimization/PowellMinimizer.cs
index a1377516..1927e44a 100644
--- a/src/Numerics/Optimization/PowellMinimizer.cs
+++ b/src/Numerics/Optimization/PowellMinimizer.cs
@@ -5,8 +5,171 @@ using System.Text;
namespace MathNet.Numerics.Optimization
{
- public class PowellSolver
+ ///
+ /// Minimizes f(p) where p is a vector of model parameters using the Powell method.
+ ///
+ public class PowellMinimizer
{
+ public double PointTolerance { get; set; }
+ public double FunctionTolerance { get; set; }
+ int? maxIterations, maxFunctionCalls;
+ double[] minimumPoint;
+ double functionAtMinimum;
+ public int FunctionCalls { get; set; }
+ public int Iterations { get; set; }
+ public int? MaxIterations { get; set; }
+ public int? MaxFunctionCalls { get; set; }
+ public double[] MinimumPoint { get { return minimumPoint; } }
+ public double FunctionValueAtMinimum { get { return functionAtMinimum; } }
+ Func function;
+ MultiDimensionalBrent powellLineSearch;
+
+ public PowellMinimizer()
+ {
+ //this.function = function;
+ //powellLineSearch = new MultiDimensionalBrent(function);
+ PointTolerance = 1e-4;
+ FunctionTolerance = 1e-4;
+ MaxIterations = null;
+ MaxFunctionCalls = null;
+ FunctionCalls = 0;
+ Iterations = 0;
+ minimumPoint = null;
+ functionAtMinimum = 0;
+ }
+
+ public double[] CurveFit(double[] x, double[] y, Func f,
+ double[] pStart)
+ {
+ // Need to minimize sum of squares of residuals; create this function:
+ Func function = (p) =>
+ {
+ double sum = 0;
+ for (int i = 0; i < x.Length; ++i)
+ {
+ double temp = y[i] - f(x[i], p);
+ sum += temp * temp;
+ }
+ return sum;
+ };
+ this.function = function;
+ powellLineSearch = new MultiDimensionalBrent(function);
+ this.Minimize(pStart);
+ return minimumPoint;
+ }
+
+ public int Minimize(double[] p)
+ {
+ // Set line search valuer to use main valuer:
+ // (this valuer takes starting point, direciton and length and
+ // returns scalar):
+ //
+ int N = p.Length; // number of dimensions
+ powellLineSearch.SetDimension(N);
+ double fval;
+ FunctionCalls = 0;
+ Iterations = 0;
+ if (maxIterations == null) maxIterations = N * 1000;
+ if (maxFunctionCalls == null) maxFunctionCalls = N * 1000;
+
+ // An array of N directions:
+ double[][] direc = new double[N][];
+ for (int i = 0; i < N; ++i)
+ {
+ direc[i] = new double[N];
+ direc[i][i] = 1.0;
+ }
+ double[] x = p;
+ double[] x1 = (double[])x.Clone();
+ powellLineSearch.Tolerance = PointTolerance * 100; // Set tolerance
+
+ fval = function(x);
+ FunctionCalls++;
+ double[] x2 = new double[N];
+ double fx;
+ double[] direc1 = new double[N];
+ double[] xnew;
+ while (true)
+ {;
+ fx = fval;
+ int bigind = 0;
+ double delta = 0.0;
+ double fx2;
+ for (int i = 0; i < N; ++i)
+ {
+ direc1 = direc[i];
+ fx2 = fval;
+ powellLineSearch.StartingPoint = x;
+ powellLineSearch.Direction = direc1;
+
+ powellLineSearch.Search(out xnew, out fval);
+ // Do a linesearch with specified starting point and direction.
+ FunctionCalls += powellLineSearch.FunctionCalls;
+
+ for (int j = 0; j < N; ++j) x[j] = xnew[j];
+
+ if ((fx2 - fval) > delta)
+ {
+ delta = fx2 - fval;
+ bigind = i;
+ }
+ }
+ Iterations++;
+ if (2.0 * (fx - fval) <= FunctionTolerance * ((Math.Abs(fx) + Math.Abs(fval)) + 1e-20)) break;
+ if (FunctionCalls >= maxFunctionCalls) break;
+ if (Iterations >= maxIterations) break;
+
+ // Construct the extrapolated point
+ direc1 = new double[N];
+ for (int i = 0; i < N; ++i)
+ {
+ direc1[i] = x[i] - x1[i];
+ x2[i] = 2.0 * x[i] - x1[i];
+ x1[i] = x[i];
+ }
+
+ fx2 = function(x2);
+ FunctionCalls++;
+ if (fx > fx2)
+ {
+ double t = 2.0 * (fx + fx2 - 2.0 * fval);
+ double temp = (fx - fval - delta);
+ t *= temp * temp;
+ temp = fx - fx2;
+ t -= delta * temp * temp;
+ if (t < 0.0)
+ {
+ powellLineSearch.StartingPoint = x;
+ powellLineSearch.Direction = direc1;
+ powellLineSearch.Search(out xnew, out fval);
+
+ FunctionCalls += powellLineSearch.FunctionCalls;
+ direc1 = new double[N];
+ for (int i = 0; i < N; ++i)
+ {
+ direc1[i] = xnew[i] - x[i];
+ x[i] = xnew[i];
+ }
+
+ direc[bigind] = direc[N - 1];
+ direc[N - 1] = direc1;
+ }
+ }
+ }
+ minimumPoint = (double[])x.Clone();
+ functionAtMinimum = fx;
+
+ // Find out what happened:
+ if (FunctionCalls >= maxFunctionCalls)
+ {
+ return 1; // Max function calls exceeded
+ }
+ if (Iterations >= maxIterations)
+ {
+ return 2; // Max iterations exceeded
+ }
+ return 0; // all good
+ }
}
}
diff --git a/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs b/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs
new file mode 100644
index 00000000..52975f64
--- /dev/null
+++ b/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs
@@ -0,0 +1,61 @@
+//
+// 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-2013 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 System;
+using MathNet.Numerics.Optimization;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+ [TestFixture]
+ public class FunctionMinimizationTests
+ {
+ [Test]
+ public void PowellCurveFit()
+ {
+ // y = b1*(1-exp[-b2*x]) + e
+ var xin = new double[] { 1, 2, 3, 5, 7, 10 };
+ var yin = new double[] { 109, 149, 149, 191, 213, 224 };
+
+ Func function = (x, p) => p[0] * (1 - Math.Exp(-p[1] * x));
+
+ var minimizer = new PowellMinimizer();
+
+ var popt = minimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }); // 100, 0.75
+
+ double[] expected = new double[] { 2.1380940889E+02, 5.4723748542E-01 };
+
+ double residual = 0;
+ for (int i = 0; i < yin.Length; ++i) residual += (yin[i] - function(xin[i], popt)) * (yin[i] - function(xin[i], popt));
+ //Assert.AreEqual(3, Brent.FindRoot(f2, 2.1, 3.4, 0.001, 50), 0.001);
+ }
+
+ }
+}
diff --git a/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs b/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs
index 325275ef..98bc6276 100644
--- a/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs
+++ b/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs
@@ -61,6 +61,5 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
for (int i = 0; i < yin.Length; ++i) residual += (yin[i] - function(xin[i], popt)) * (yin[i] - function(xin[i], popt));
//Assert.AreEqual(3, Brent.FindRoot(f2, 2.1, 3.4, 0.001, 50), 0.001);
}
-
}
}
diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj
index 0776c7da..8c2292b6 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -353,6 +353,7 @@
+