Cleaning up Powell and Brent minimizers.

src/Numerics/Optimization/BrentMinimizer.cs

@@ -11,167 +11,80 @@ namespace MathNet.Numerics.Optimization
     public class BrentOptions
     {
         public int MaximumIterations = 1000;
-
         public double FunctionTolerance = 1e-4;
     }
+
+    /// <summary>
+    /// Result of Brent Minimization.
+    /// </summary>
+    public class BrentResult
+    {
+        public int NumberOfIterations;
+        public double MinimumPoint;
+        public double MinimumFunctionValue;
+    }
     
     /// <summary>
     /// Minimizes f(p) where p is a model parameter scalar, i.e. a line-search.
+    /// Inspired by the SciPy implementation.
     /// </summary>
     public class BrentMinimizer
     {
-        const double verySmallNumber = 1e-21, goldenRatio = 1.618034, minimumTolerance = 1.0e-11;
+        public struct Bracket
-        const double growLimit = 110.0, conjugateGradient = 0.3819660;
-        const int maxIterations = 500;
-
-        int bracketFunctionCalls;
-        public int FunctionCalls, Iterations;
-        double minPoint, minFunction;
-        Func<double, double> function;
-        double pointA, pointB, pointC;
-        double functionA, functionB, functionC;
-
-        public double Tolerance { get; set; }
-
-        public BrentMinimizer(Func<double, double> function)
         {
-            this.function = function;
+            public double PointA;
-            Tolerance = 1e-4;
+            public double PointB;
+            public double PointC;
+            public double FunctionA;
+            public double FunctionB;
+            public double FunctionC;
         }
+        
+        public BrentResult Result { get; private set; }
 
-        public int Search(out double minPoint, out double minFunction)
+        public readonly BrentOptions Options = new BrentOptions();
-        {
+        
-            bracketFunctionCalls = 0;
+        const double verySmallNumber = 1e-21, goldenRatio = 1.618034, minimumTolerance = 1.0e-11;
-            UpdateBracketInterval(0, 1);
+        const double growLimit = 110.0, conjugateGradient = 0.3819660;
-            FunctionCalls = 0;
-            Iterations = 0;
-            int result = BrentMinimize();
-            FunctionCalls += bracketFunctionCalls;
-            minPoint = this.minPoint;
-            minFunction = this.minFunction;
-            return result;
-        }
 
-        private int UpdateBracketInterval(double pointAStart, double pointBStart)
+        /// <summary>
+        /// Find the minimum of the supplied function using the Brent method.
+        /// </summary>
+        /// <param name="function"></param>
+        /// <returns></returns>
+        public double Minimize(Func<double, double> function)
         {
+            Bracket bracket;
+            UpdateBracketInterval(function, new Bracket() { PointA = 0, PointB = 1 }, out bracket);
             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
+            x = w = v = bracket.PointB; // x is the point with lowest function value encountered
             fw = fv = fx = function(x);
-            if (pointA < pointC)
+            if (bracket.PointA < bracket.PointC)
             {
-                a = pointA; b = pointC;
+                a = bracket.PointA; b = bracket.PointC;
             }
             else
             {
-                a = pointC; b = pointA;
+                a = bracket.PointC; b = bracket.PointA;
             }
             deltax = 0.0;
-            FunctionCalls = 1;
-            Iterations = 0;
             rat = 0;
-            while (Iterations < maxIterations)
+            while (iterations < Options.MaximumIterations)
             {
-                double tol1, tol2, xmin, fval, xmid;
+                double tol1, tol2, xmid;
                 double temp1, temp2, p;
                 double u, fu, dx_temp;
-                tol1 = Tolerance * Math.Abs(x) + minimumTolerance;
+                tol1 = Options.FunctionTolerance * 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 
+                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 
+                    if (x >= xmid) deltax = a - x; // do a golden section step 
                     else deltax = b - x;
                     rat = cg * deltax;
                 }
@@ -215,7 +128,6 @@ namespace MathNet.Numerics.Optimization
                     u = x + rat;
                 }
                 fu = function(u);
-                FunctionCalls++;
                 if (fu > fx)        // if it's bigger than current 
                 {
                     if (u < x) a = u;
@@ -236,70 +148,127 @@ namespace MathNet.Numerics.Optimization
                     v = w; w = x; x = u;
                     fv = fw; fw = fx; fx = fu;
                 }
-                Iterations++;
+                iterations++;
             }
-            this.minPoint = x;
+            this.Result = new BrentResult() { NumberOfIterations = iterations, MinimumPoint = x, MinimumFunctionValue = fx };
-            this.minFunction = fx;
+            return x;            
-            return result;
-        }
-    }
-
-    /// <summary>
-    /// Minimizes f(p u) where p is a model parameter scalar and u is a direction vector.
-    /// </summary>
-    public class MultiDimensionalBrent
-    {
-        private Func<double[], double> function;
-        private int functionCalls;
-        double[] point;
-        BrentMinimizer lineSearch;
-
-        public MultiDimensionalBrent(Func<double[], double> 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
+        /// <summary>
+        /// Find the minimum of the supplied function along a specified line, using the Brent method.
+        /// </summary>
+        /// <param name="function"></param>
+        /// <param name="direction">Direction of line.</param>
+        /// <param name="startingPoint">Starting point of line.</param>
+        /// <returns></returns>
+        public double Minimize(Func<double[], double> function, double[] direction, double[] startingPoint, out double[] minimumPoint)
         {
-            get { return lineSearch.FunctionCalls; }
+            double[] point = new double[direction.Length];
-        }
+            Func<double, double> functionAlongLine = (p) =>
-
+                {
-        public double[] StartingPoint { get; set; }
+                    for (int i = 0; i < point.Length; ++i)
-        public double[] Direction { get; set; }
+                        point[i] = startingPoint[i] + direction[i] * p;
-
+                    return function(point);
-
+                };
-        public void SetDimension(int N)
+            double result = Minimize(functionAlongLine);
-        {
+            minimumPoint = point;
-            point = new double[N];
+            return result;
         }
 
-        public int Search(out double[] minPoint, out double minFunction)
+        /// <summary>
+        /// Updates the bracket.
+        /// </summary>
+        /// <param name="function"></param>
+        /// <param name="bracketInitial"></param>
+        /// <param name="newBracket"></param>
+        /// <returns></returns>
+        private static bool UpdateBracketInterval(Func<double, double> function, Bracket bracketInitial, out Bracket newBracket)
         {
-            double path;
+            int iterations = 0;
-            int result = lineSearch.Search(out path, out minFunction);
+            double pointA = bracketInitial.PointA;
-            for (int i = 0; i < StartingPoint.Length; ++i)
+            double pointB = bracketInitial.PointB;
+            int maxIterations;
+            maxIterations = 1000;
+            double functionA = function(pointA);
+            double functionB = function(pointB);
+            double temp;
+            if (functionA < functionB) // Swap points over 
             {
-                point[i] = StartingPoint[i] + path * Direction[i];
+                temp = functionA; functionA = functionB; functionB = temp;
+                temp = pointA; pointA = pointB; pointB = temp;
             }
-            minPoint = point;
+            double pointC = pointB + goldenRatio * (pointB - pointA);
-            return result;
+            double functionC = function(pointC);
-        }
+            iterations = 0;
-
+            double temp1, temp2, value, denom;
-        public double PointAlongLine(double path)
+            while (functionC < functionB)
-        {
-            for (int i = 0; i < StartingPoint.Length; ++i)
             {
-                point[i] = StartingPoint[i] + path * Direction[i];
+                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)
+                {
+                    newBracket = bracketInitial;
+                    return false;
+                }
+                iterations++;
+                if ((pointW - pointC) * (pointB - pointW) > 0.0)
+                {
+                    functionW = function(pointW);
+                    if (functionW < functionC)
+                    {
+                        pointA = pointB; pointB = pointW;
+                        functionA = functionB; functionB = functionW;
+                        break;
+                    }
+                    else if (functionW > functionB)
+                    {
+                        pointC = pointW; functionC = functionW;
+                        break;
+                    }
+                    pointW = pointC + goldenRatio * (pointC - pointB);
+                    functionW = function(pointW);
+                }
+                else if ((pointW - wlim) * (wlim - pointC) >= 0.0)
+                {
+                    pointW = wlim;
+                    functionW = function(pointW);
+                }
+                else if ((pointW - wlim) * (pointC - pointW) > 0.0)
+                {
+                    functionW = function(pointW);
+                    if (functionW < functionC)
+                    {
+                        pointB = pointC; pointC = pointW;
+                        pointW = pointC + goldenRatio * (pointC - pointB);
+                        functionB = functionC; functionC = functionW;
+                        functionW = function(pointW);
+                    }
+                }
+                else
+                {
+                    pointW = pointC + goldenRatio * (pointC - pointB);
+                    functionW = function(pointW);
+                }
+                pointA = pointB; pointB = pointC; pointC = pointW;
+                functionA = functionB; functionB = functionC; functionC = functionW;
             }
-            return function(point);
+            newBracket = new Bracket()
+            {
+                PointA = pointA,
+                PointB = pointB,
+                PointC = pointC,
+                FunctionA = functionA,
+                FunctionB = functionB,
+                FunctionC = functionC
+            };
+            return true;
         }
+
     }
 }

src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs

@@ -13,9 +13,7 @@ namespace MathNet.Numerics.Optimization
     public class NonLinearLeastSquaresOptions
     {
         public int MaximumIterations = 1000;
-
         public int MaximumTrialStepIterations = 100;
-
         public NonLinearLeastSquaresConvergenceType ConvergenceType;
 
         /// <summary>
@@ -54,7 +52,7 @@ namespace MathNet.Numerics.Optimization
     /// <summary>
     /// For details of convergence criteria, see Options.
     /// </summary>
-    public enum NonLinearLeastSquaresConvergenceType { NoneMaxIterationExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, Error };
+    public enum NonLinearLeastSquaresConvergenceType { MaxIterationsExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, Error };
 
     /// <summary>
     /// Result of Non-Linear Least Squares Minimization.
@@ -62,7 +60,6 @@ namespace MathNet.Numerics.Optimization
     public class NonLinearLeastSquaresResult
     {
         public int NumberOfIterations;
-
         public NonLinearLeastSquaresConvergenceType ConvergenceType;
     }
 

src/Numerics/Optimization/PowellMinimizer.cs

@@ -5,39 +5,39 @@ using System.Text;
 
 namespace MathNet.Numerics.Optimization
 {
+    /// <summary>
+    /// Options for Powell Minimization.
+    /// </summary>
+    public class PowellOptions
+    {
+        public int? MaximumIterations = null;
+        public int? MaximumFunctionCalls = null;
+        public double PointTolerance = 1e-4;
+        public double FunctionTolerance = 1e-4;
+    }
+
+    public enum PowellConvergenceType { Success, MaxIterationsExceeded, MaxFunctionCallsExceeded };
+
+    /// <summary>
+    /// Result of Powell Minimization.
+    /// </summary>
+    public class PowellResult
+    {
+        public int NumberOfIterations;
+        public int NumberOfFunctionCalls;
+        public double[] MinimumPoint;
+        public double MinimumFunctionValue;
+        public PowellConvergenceType ConvergenceType; 
+    }
+    
     /// <summary>
     /// Minimizes f(p) where p is a vector of model parameters using the Powell method.
     /// </summary>
     public class PowellMinimizer
     {
-        public double PointTolerance { get; set; }
+        public PowellResult Result { get; private set; }
-        public double FunctionTolerance { get; set; }
+
-        int? maxIterations, maxFunctionCalls;
+        public readonly PowellOptions Options = new PowellOptions();
-        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<double[], double> 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<double, double[], double> f,
             double[] pStart)
@@ -53,61 +53,71 @@ namespace MathNet.Numerics.Optimization
                 }
                 return sum;
             };
-            this.function = function;
+            return Minimize(function, pStart);
-            powellLineSearch = new MultiDimensionalBrent(function);
-            this.Minimize(pStart);
-            return minimumPoint;
         }
 
-        public int Minimize(double[] p)
+        public double[] Minimize(Func<double[], double> function, double[] p)
         {
-            // Set line search valuer to use main valuer: 
+            BrentMinimizer brentMinimizer = new BrentMinimizer();
-            // (this valuer takes starting point, direciton and length and 
+            // used in closure:
-            // returns scalar): 
+            double[] point = new double[p.Length];
-            // 
+            double[] startingPoint = new double[p.Length];
-            int N = p.Length; // number of dimensions 
+            double[] direction = new double[p.Length];
-            powellLineSearch.SetDimension(N);
+            double lineMiniumum = 0;
+            int functionCalls = 0;
+            Func<double, double> functionAlongLine = (u) =>
+            {
+                for (int i = 0; i < point.Length; ++i)
+                    point[i] = startingPoint[i] + direction[i] * u;
+                lineMiniumum = function(point);
+                functionCalls++;
+                return lineMiniumum;
+            };
+
+            int n = p.Length; // number of dimensions 
             double fval;
-            FunctionCalls = 0;
+
-            Iterations = 0;
+            int iterations = 0;
-            if (maxIterations == null) maxIterations = N * 1000;
+            int maxIterations = (Options.MaximumIterations == null) ? n * 1000 : (int)Options.MaximumIterations;
-            if (maxFunctionCalls == null) maxFunctionCalls = N * 1000;
+            int maxFunctionCalls = (Options.MaximumFunctionCalls == null) ? n * 1000 : (int)Options.MaximumFunctionCalls;
-
+
-            // An array of N directions:
+            // An array of n directions:
-            double[][] direc = new double[N][];
+            double[][] direc = new double[n][];
-            for (int i = 0; i < N; ++i)
+            for (int i = 0; i < n; ++i)
             {
-                direc[i] = new double[N];
+                direc[i] = new double[n];
                 direc[i][i] = 1.0;
             }
             double[] x = p;
             double[] x1 = (double[])x.Clone(); 
-            powellLineSearch.Tolerance = PointTolerance * 100; // Set tolerance 
+            
+            brentMinimizer.Options.FunctionTolerance = Options.PointTolerance * 100;
 
             fval = function(x);
-            FunctionCalls++;
+
-            double[] x2 = new double[N];
+            double[] x2 = new double[n];
             double fx;
-            double[] direc1 = new double[N];
+            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)
+                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.
-                    // Do a linesearch with specified starting point and direction. 
+                    direction = direc1;
-                    FunctionCalls += powellLineSearch.FunctionCalls;
+                    startingPoint = x;
+                    double u = brentMinimizer.Minimize(functionAlongLine);
+                    fval = functionAlongLine(u);
+                    xnew = point; 
 
-                    for (int j = 0; j < N; ++j) x[j] = xnew[j];
+                    for (int j = 0; j < n; ++j) x[j] = xnew[j];
 
                     if ((fx2 - fval) > delta)
                     {
@@ -115,14 +125,14 @@ namespace MathNet.Numerics.Optimization
                         bigind = i;
                     }
                 }
-                Iterations++;
+                iterations++;
-                if (2.0 * (fx - fval) <= FunctionTolerance * ((Math.Abs(fx) + Math.Abs(fval)) + 1e-20)) break;
+                if (2.0 * (fx - fval) <= Options.FunctionTolerance * ((Math.Abs(fx) + Math.Abs(fval)) + 1e-20)) break;
-                if (FunctionCalls >= maxFunctionCalls) break;
+                if (functionCalls >= maxFunctionCalls) break;
-                if (Iterations >= maxIterations) break;
+                if (iterations >= maxIterations) break;
 
                 // Construct the extrapolated point  
-                direc1 = new double[N];
+                direc1 = new double[n];
-                for (int i = 0; i < N; ++i)
+                for (int i = 0; i < n; ++i)
                 {
                     direc1[i] = x[i] - x1[i];
                     x2[i] = 2.0 * x[i] - x1[i];
@@ -130,7 +140,6 @@ namespace MathNet.Numerics.Optimization
                 }
 
                 fx2 = function(x2);
-                FunctionCalls++;
                 if (fx > fx2)
                 {
                     double t = 2.0 * (fx + fx2 - 2.0 * fval);
@@ -140,36 +149,40 @@ namespace MathNet.Numerics.Optimization
                     t -= delta * temp * temp;
                     if (t < 0.0)
                     {
-                        powellLineSearch.StartingPoint = x;
+                        direction = direc1;
-                        powellLineSearch.Direction = direc1;
+                        startingPoint = x;
-                        powellLineSearch.Search(out xnew, out fval);
+                        double u = brentMinimizer.Minimize(functionAlongLine);
+                        fval = functionAlongLine(u);
+                        xnew = point;
 
-                        FunctionCalls += powellLineSearch.FunctionCalls;
+                        direc1 = new double[n];
-                        direc1 = new double[N];
+                        for (int i = 0; i < n; ++i)
-                        for (int i = 0; i < N; ++i)
                         {
                             direc1[i] = xnew[i] - x[i];
                             x[i] = xnew[i];
                         }
 
-                        direc[bigind] = direc[N - 1];
+                        direc[bigind] = direc[n - 1];
-                        direc[N - 1] = direc1;
+                        direc[n - 1] = direc1;
                     }
                 }
             }
-            minimumPoint = (double[])x.Clone();
+            var convergenceType = PowellConvergenceType.Success;
-            functionAtMinimum = fx;
+            if (functionCalls >= maxFunctionCalls) 
-            
+                convergenceType = PowellConvergenceType.MaxFunctionCallsExceeded;
-            // Find out what happened: 
+            else if (iterations > maxIterations) 
-            if (FunctionCalls >= maxFunctionCalls)
+                convergenceType = PowellConvergenceType.MaxFunctionCallsExceeded;
-            {
+
-                return 1; // Max function calls exceeded 
+            Result = new PowellResult()
-            }
-            if (Iterations >= maxIterations)
             {
-                return 2; // Max iterations exceeded 
+                MinimumPoint = (double[])x.Clone(),
-            }
+                MinimumFunctionValue = fx,
-            return 0; // all good 
+                ConvergenceType = convergenceType,
+                NumberOfIterations = iterations,
+                NumberOfFunctionCalls = functionCalls
+            };
+            
+            return Result.MinimumPoint;
         }
     }
 }

src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs

@@ -180,7 +180,7 @@ namespace MathNet.Numerics.Providers.Optimization.Mkl
             switch (rciRequest)
             {
                 case -1:
-                    convergenceType = NonLinearLeastSquaresConvergenceType.NoneMaxIterationExceeded; break;
+                    convergenceType = NonLinearLeastSquaresConvergenceType.MaxIterationsExceeded; break;
                 case -2:
                     convergenceType = NonLinearLeastSquaresConvergenceType.Criterion0; break;
                 case -3:

src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs

@@ -48,14 +48,22 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
 
             var minimizer = new PowellMinimizer();
 
-            var popt = minimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }); // 100, 0.75
+            var watch = new System.Diagnostics.Stopwatch(); watch.Start();
+            double[] popt = null;
+            for (int i = 0; i < 1000; ++i)
+            {
+                popt = minimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }); // 100, 0.75
+            }
+            watch.Stop();
+            double elapsed = watch.ElapsedMilliseconds;
 
             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);
-        }
 
+            Assert.AreEqual(expected[0], popt[0], 1e-4);
+            Assert.AreEqual(expected[1], popt[1], 1e-4);
+        }
     }
 }

src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs

@@ -46,20 +46,24 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
             // 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 };
-            var popt = minimizer.CurveFit(xin, yin, (x, p) => p[0] * (1 - Math.Exp(-p[1] * x)), new double[] { 1, 1 });
+            
+            // estimated derivative method: does not find best solution.
+            //var popt = minimizer.CurveFit(xin, yin, (x, p) => p[0] * (1 - Math.Exp(-p[1] * x)), new double[] { 1, 1 });
 
             Func<double, double[], double> function = (x, p) => p[0] * (1 - Math.Exp(-p[1] * x));
             Func<double, double[], double[]> jacobian = (x, p) => new double[] { 
                 1 - Math.Exp(-p[1] * x), 
                 p[0] * x * Math.Exp(-p[1] * x) };
 
-            popt = minimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }, jacobian); // 100, 0.75
+            var popt = minimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }, jacobian); // 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)); 
+            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);
+
+            Assert.AreEqual(expected[0], popt[0], 1e-6);
+            Assert.AreEqual(expected[1], popt[1], 1e-6);
         }
     }
 }