diff --git a/src/Numerics/Optimization/BrentMinimizer.cs b/src/Numerics/Optimization/BrentMinimizer.cs index 3d029903..78d3c6f7 100644 --- a/src/Numerics/Optimization/BrentMinimizer.cs +++ b/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; } + + /// + /// Result of Brent Minimization. + /// + public class BrentResult + { + public int NumberOfIterations; + public double MinimumPoint; + public double MinimumFunctionValue; + } /// /// Minimizes f(p) where p is a model parameter scalar, i.e. a line-search. + /// Inspired by the SciPy implementation. /// 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) + public struct Bracket { - this.function = function; - Tolerance = 1e-4; + public double PointA; + 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) - { - bracketFunctionCalls = 0; - UpdateBracketInterval(0, 1); - FunctionCalls = 0; - Iterations = 0; - int result = BrentMinimize(); - FunctionCalls += bracketFunctionCalls; - minPoint = this.minPoint; - minFunction = this.minFunction; - return result; - } + public readonly BrentOptions Options = new BrentOptions(); + + const double verySmallNumber = 1e-21, goldenRatio = 1.618034, minimumTolerance = 1.0e-11; + const double growLimit = 110.0, conjugateGradient = 0.3819660; - private int UpdateBracketInterval(double pointAStart, double pointBStart) + /// + /// Find the minimum of the supplied function using the Brent method. + /// + /// + /// + public double Minimize(Func 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 - { - xmin = x; fval = fx; - result = 1; + if (Math.Abs(x - xmid) < (tol2 - 0.5 * (b - a))) // check for convergence 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.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; } + this.Result = new BrentResult() { NumberOfIterations = iterations, MinimumPoint = x, MinimumFunctionValue = fx }; + return x; } - public int FunctionCalls + /// + /// Find the minimum of the supplied function along a specified line, using the Brent method. + /// + /// + /// Direction of line. + /// Starting point of line. + /// + public double Minimize(Func function, double[] direction, double[] startingPoint, out double[] minimumPoint) { - get { return lineSearch.FunctionCalls; } - } - - public double[] StartingPoint { get; set; } - public double[] Direction { get; set; } - - - public void SetDimension(int N) - { - point = new double[N]; + double[] point = new double[direction.Length]; + Func functionAlongLine = (p) => + { + for (int i = 0; i < point.Length; ++i) + point[i] = startingPoint[i] + direction[i] * p; + return function(point); + }; + double result = Minimize(functionAlongLine); + minimumPoint = point; + return result; } - public int Search(out double[] minPoint, out double minFunction) + /// + /// Updates the bracket. + /// + /// + /// + /// + /// + private static bool UpdateBracketInterval(Func function, Bracket bracketInitial, out Bracket newBracket) { - double path; - int result = lineSearch.Search(out path, out minFunction); - for (int i = 0; i < StartingPoint.Length; ++i) + int iterations = 0; + double pointA = bracketInitial.PointA; + 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; - return result; - } - - public double PointAlongLine(double path) - { - for (int i = 0; i < StartingPoint.Length; ++i) + double pointC = pointB + goldenRatio * (pointB - pointA); + double functionC = function(pointC); + iterations = 0; + double temp1, temp2, value, denom; + while (functionC < functionB) { - 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; } + } } diff --git a/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs b/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs index bec4f902..caaf1a99 100644 --- a/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs +++ b/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; /// @@ -54,7 +52,7 @@ namespace MathNet.Numerics.Optimization /// /// For details of convergence criteria, see Options. /// - public enum NonLinearLeastSquaresConvergenceType { NoneMaxIterationExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, Error }; + public enum NonLinearLeastSquaresConvergenceType { MaxIterationsExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, Error }; /// /// Result of Non-Linear Least Squares Minimization. @@ -62,7 +60,6 @@ namespace MathNet.Numerics.Optimization public class NonLinearLeastSquaresResult { public int NumberOfIterations; - public NonLinearLeastSquaresConvergenceType ConvergenceType; } diff --git a/src/Numerics/Optimization/PowellMinimizer.cs b/src/Numerics/Optimization/PowellMinimizer.cs index 1927e44a..2d062d04 100644 --- a/src/Numerics/Optimization/PowellMinimizer.cs +++ b/src/Numerics/Optimization/PowellMinimizer.cs @@ -5,39 +5,39 @@ using System.Text; namespace MathNet.Numerics.Optimization { + /// + /// Options for Powell Minimization. + /// + 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 }; + + /// + /// Result of Powell Minimization. + /// + public class PowellResult + { + public int NumberOfIterations; + public int NumberOfFunctionCalls; + public double[] MinimumPoint; + public double MinimumFunctionValue; + public PowellConvergenceType ConvergenceType; + } + /// /// 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 PowellResult Result { get; private set; } + + public readonly PowellOptions Options = new PowellOptions(); public double[] CurveFit(double[] x, double[] y, Func f, double[] pStart) @@ -53,61 +53,71 @@ namespace MathNet.Numerics.Optimization } return sum; }; - this.function = function; - powellLineSearch = new MultiDimensionalBrent(function); - this.Minimize(pStart); - return minimumPoint; + return Minimize(function, pStart); } - public int Minimize(double[] p) + public double[] Minimize(Func function, 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); + BrentMinimizer brentMinimizer = new BrentMinimizer(); + // used in closure: + double[] point = new double[p.Length]; + double[] startingPoint = new double[p.Length]; + double[] direction = new double[p.Length]; + double lineMiniumum = 0; + int functionCalls = 0; + Func 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; - 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) + + int iterations = 0; + int maxIterations = (Options.MaximumIterations == null) ? n * 1000 : (int)Options.MaximumIterations; + int maxFunctionCalls = (Options.MaximumFunctionCalls == null) ? n * 1000 : (int)Options.MaximumFunctionCalls; + + // An array of n directions: + double[][] direc = new double[n][]; + 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. - FunctionCalls += powellLineSearch.FunctionCalls; + // Do a linesearch with specified starting point and direction. + direction = direc1; + 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++; - if (2.0 * (fx - fval) <= FunctionTolerance * ((Math.Abs(fx) + Math.Abs(fval)) + 1e-20)) break; - if (FunctionCalls >= maxFunctionCalls) break; - if (Iterations >= maxIterations) break; + iterations++; + if (2.0 * (fx - fval) <= Options.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 = new double[n]; + 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; - powellLineSearch.Direction = direc1; - powellLineSearch.Search(out xnew, out fval); + direction = direc1; + startingPoint = x; + double u = brentMinimizer.Minimize(functionAlongLine); + fval = functionAlongLine(u); + xnew = point; - FunctionCalls += powellLineSearch.FunctionCalls; - direc1 = new double[N]; - for (int i = 0; i < N; ++i) + 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; + 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) + var convergenceType = PowellConvergenceType.Success; + if (functionCalls >= maxFunctionCalls) + convergenceType = PowellConvergenceType.MaxFunctionCallsExceeded; + else if (iterations > maxIterations) + convergenceType = PowellConvergenceType.MaxFunctionCallsExceeded; + + Result = new PowellResult() { - return 2; // Max iterations exceeded - } - return 0; // all good + MinimumPoint = (double[])x.Clone(), + MinimumFunctionValue = fx, + ConvergenceType = convergenceType, + NumberOfIterations = iterations, + NumberOfFunctionCalls = functionCalls + }; + + return Result.MinimumPoint; } } } diff --git a/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs b/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs index b5530720..d9f81a63 100644 --- a/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs +++ b/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: diff --git a/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs b/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs index 52975f64..14b5cb16 100644 --- a/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs +++ b/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); + } } } diff --git a/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs b/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs index 98bc6276..5c345a44 100644 --- a/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs +++ b/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 function = (x, p) => p[0] * (1 - Math.Exp(-p[1] * x)); Func 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)); - //Assert.AreEqual(3, Brent.FindRoot(f2, 2.1, 3.4, 0.001, 50), 0.001); + for (int i = 0; i < yin.Length; ++i) residual += (yin[i] - function(xin[i], popt)) * (yin[i] - function(xin[i], popt)); + + Assert.AreEqual(expected[0], popt[0], 1e-6); + Assert.AreEqual(expected[1], popt[1], 1e-6); } } }