diff --git a/src/Numerics/Optimization/BrentMinimizer.cs b/src/Numerics/Optimization/BrentMinimizer.cs index 8cd5eb7b..fab72481 100644 --- a/src/Numerics/Optimization/BrentMinimizer.cs +++ b/src/Numerics/Optimization/BrentMinimizer.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -50,7 +50,7 @@ namespace MathNet.Numerics.Optimization 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. @@ -66,13 +66,13 @@ namespace MathNet.Numerics.Optimization public double FunctionB; public double FunctionC; } - + public BrentResult Result { get; private set; } 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; + + const double VerySmallNumber = 1e-21, GoldenRatio = 1.618034, MinimumTolerance = 1.0e-11; + const double GrowLimit = 110.0, ConjugateGradient = 0.3819660; /// /// Find the minimum of the supplied function using the Brent method. @@ -82,20 +82,22 @@ namespace MathNet.Numerics.Optimization public double Minimize(Func function) { Bracket bracket; - UpdateBracketInterval(function, new Bracket() { PointA = 0, PointB = 1 }, out bracket); + UpdateBracketInterval(function, new Bracket { PointA = 0, PointB = 1 }, out bracket); int iterations = 0; double x, w, v, fx, fw, fv; double a, b, deltax, rat; - double cg = conjugateGradient; + double cg = ConjugateGradient; x = w = v = bracket.PointB; // x is the point with lowest function value encountered fw = fv = fx = function(x); if (bracket.PointA < bracket.PointC) { - a = bracket.PointA; b = bracket.PointC; + a = bracket.PointA; + b = bracket.PointC; } else { - a = bracket.PointC; b = bracket.PointA; + a = bracket.PointC; + b = bracket.PointA; } deltax = 0.0; rat = 0; @@ -104,32 +106,32 @@ namespace MathNet.Numerics.Optimization double tol1, tol2, xmid; double temp1, temp2, p; double u, fu, dx_temp; - tol1 = Options.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 + tol1 = Options.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 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; + rat = cg*deltax; } - else // do a parabolic step + 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); + 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))) + // 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. + rat = p*1.0/temp2; // if parabolic step is useful. u = x + rat; if (((u - a) < tol2) || ((b - u) < tol2)) { @@ -139,13 +141,13 @@ namespace MathNet.Numerics.Optimization } else { - if (x >= xmid) deltax = a - x; // if it's not do a golden section step + if (x >= xmid) deltax = a - x; // if it's not do a golden section step else deltax = b - x; - rat = cg * deltax; + rat = cg*deltax; } } - if (Math.Abs(rat) < tol1) // update by at least tol1 + if (Math.Abs(rat) < tol1) // update by at least tol1 { if (rat >= 0) u = x + tol1; else u = x - tol1; @@ -155,30 +157,38 @@ namespace MathNet.Numerics.Optimization u = x + rat; } fu = function(u); - if (fu > fx) // if it's bigger than current + 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; + v = w; + w = u; + fv = fw; + fw = fu; } else if ((fu <= fv) || (v == x) || (v == w)) { - v = u; fv = fu; + 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; + v = w; + w = x; + x = u; + fv = fw; + fw = fx; + fx = fu; } iterations++; } - this.Result = new BrentResult() { NumberOfIterations = iterations, MinimumPoint = x, MinimumFunctionValue = fx }; - return x; + Result = new BrentResult { NumberOfIterations = iterations, MinimumPoint = x, MinimumFunctionValue = fx }; + return x; } /// @@ -191,12 +201,12 @@ namespace MathNet.Numerics.Optimization public double Minimize(Func function, double[] direction, double[] startingPoint, out double[] minimumPoint) { 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); - }; + 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; @@ -209,83 +219,93 @@ namespace MathNet.Numerics.Optimization /// /// /// - private static bool UpdateBracketInterval(Func function, Bracket bracketInitial, out Bracket newBracket) + static bool UpdateBracketInterval(Func function, Bracket bracketInitial, out Bracket newBracket) { - int iterations = 0; double pointA = bracketInitial.PointA; double pointB = bracketInitial.PointB; - int maxIterations; - maxIterations = 1000; + const int maxIterations = 1000; double functionA = function(pointA); double functionB = function(pointB); - double temp; - if (functionA < functionB) // Swap points over + if (functionA < functionB) // Swap points over { - temp = functionA; functionA = functionB; functionB = temp; - temp = pointA; pointA = pointB; pointB = temp; + double temp = functionA; + functionA = functionB; + functionB = temp; + temp = pointA; + pointA = pointB; + pointB = temp; } - double pointC = pointB + goldenRatio * (pointB - pointA); + double pointC = pointB + GoldenRatio*(pointB - pointA); double functionC = function(pointC); - iterations = 0; + int iterations = 0; double temp1, temp2, value, denom; while (functionC < functionB) { - double pointW, functionW, wlim; - temp1 = (pointB - pointA) * (functionB - functionC); - temp2 = (pointB - pointC) * (functionB - functionA); + double functionW; + 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 (Math.Abs(value) < VerySmallNumber) denom = 2.0*VerySmallNumber; + else denom = 2.0*value; + double pointW = pointB - ((pointB - pointC)*temp2 - (pointB - pointA)*temp1)/denom; + double wlim = pointB + GrowLimit*(pointC - pointB); if (iterations > maxIterations) { newBracket = bracketInitial; return false; } iterations++; - if ((pointW - pointC) * (pointB - pointW) > 0.0) + if ((pointW - pointC)*(pointB - pointW) > 0.0) { functionW = function(pointW); if (functionW < functionC) { - pointA = pointB; pointB = pointW; - functionA = functionB; functionB = functionW; + pointA = pointB; + pointB = pointW; + functionA = functionB; + functionB = functionW; break; } - else if (functionW > functionB) + if (functionW > functionB) { - pointC = pointW; functionC = functionW; + pointC = pointW; + functionC = functionW; break; } - pointW = pointC + goldenRatio * (pointC - pointB); + pointW = pointC + GoldenRatio*(pointC - pointB); functionW = function(pointW); } - else if ((pointW - wlim) * (wlim - pointC) >= 0.0) + else if ((pointW - wlim)*(wlim - pointC) >= 0.0) { pointW = wlim; functionW = function(pointW); } - else if ((pointW - wlim) * (pointC - pointW) > 0.0) + 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; + pointB = pointC; + pointC = pointW; + pointW = pointC + GoldenRatio*(pointC - pointB); + functionB = functionC; + functionC = functionW; functionW = function(pointW); } } else { - pointW = pointC + goldenRatio * (pointC - pointB); + pointW = pointC + GoldenRatio*(pointC - pointB); functionW = function(pointW); } - pointA = pointB; pointB = pointC; pointC = pointW; - functionA = functionB; functionB = functionC; functionC = functionW; + pointA = pointB; + pointB = pointC; + pointC = pointW; + functionA = functionB; + functionB = functionC; + functionC = functionW; } - newBracket = new Bracket() + newBracket = new Bracket { PointA = pointA, PointB = pointB, diff --git a/src/Numerics/Optimization/IMinimizer.cs b/src/Numerics/Optimization/IMinimizer.cs index 3e333a11..0515468d 100644 --- a/src/Numerics/Optimization/IMinimizer.cs +++ b/src/Numerics/Optimization/IMinimizer.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -28,13 +28,13 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; + namespace MathNet.Numerics.Optimization { - using System; - /// /// Interface implemented by a class that minimizes f(p) where p is a vector of model parameters. - /// A class implenting this interface can then be used to solve curve fitting problems. + /// A class implenting this interface can then be used to solve curve fitting problems. /// public interface IMinimizer { @@ -56,13 +56,13 @@ namespace MathNet.Numerics.Optimization double[] pStart) { // Need to minimize sum of squares of residuals; create this function: - Func function = (p) => + 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; + sum += temp*temp; } return sum; }; diff --git a/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs b/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs index dbf8d739..9757790d 100644 --- a/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs +++ b/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -47,7 +47,7 @@ namespace MathNet.Numerics.Optimization public NonLinearLeastSquaresConvergenceType ConvergenceType; /// - /// Convergence if Δ < Criterion0, Δ is trust region size. + /// Convergence if Δ < Criterion0, Δ is trust region size. /// public double Criterion0 = 1e-7; @@ -82,7 +82,16 @@ namespace MathNet.Numerics.Optimization /// /// For details of convergence criteria, see Options. /// - public enum NonLinearLeastSquaresConvergenceType { MaxIterationsExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, Error }; + public enum NonLinearLeastSquaresConvergenceType + { + MaxIterationsExceeded, + Criterion0, + Criterion1, + Criterion2, + Criterion3, + Criterion4, + Error + }; /// /// Result of Non-Linear Least Squares Minimization. @@ -97,7 +106,7 @@ namespace MathNet.Numerics.Optimization /// /// This class is a special function minimizer that minimizes functions of the form - /// f(p) = |r(p)|^2 where r is a vector of residuals and p is a vector of model parameters. + /// f(p) = |r(p)|^2 where r is a vector of residuals and p is a vector of model parameters. /// public class NonLinearLeastSquaresMinimizer // Note does not implement IMinimizer, since it curve fitting problems can be solved more efficiently. { @@ -129,15 +138,16 @@ namespace MathNet.Numerics.Optimization // jac is df_i / dp_j Jacobian jacobianFunction = null; - if (jacobian != null) jacobianFunction = (p, jac) => - { - for (int i = 0; i < y.Length; ++i) + if (jacobian != null) + jacobianFunction = (p, jac) => { - double[] values = jacobian(x[i], p); - for (int j = 0; j < values.Length; ++j) - jac[j * y.Length + i] = -values[j]; - } - }; + for (int i = 0; i < y.Length; ++i) + { + double[] values = jacobian(x[i], p); + for (int j = 0; j < values.Length; ++j) + jac[j*y.Length + i] = -values[j]; + } + }; double[] parameters; Result = provider.NonLinearLeastSquaresUnboundedMinimize(y.Length, pStart, function, out parameters, jacobianFunction); diff --git a/src/Numerics/Optimization/PowellMinimizer.cs b/src/Numerics/Optimization/PowellMinimizer.cs index 4bb5931c..8be05233 100644 --- a/src/Numerics/Optimization/PowellMinimizer.cs +++ b/src/Numerics/Optimization/PowellMinimizer.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -28,11 +28,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; + namespace MathNet.Numerics.Optimization { - using System; - using System.Linq; - /// /// Options for Powell Minimization. /// @@ -44,7 +43,12 @@ namespace MathNet.Numerics.Optimization public double FunctionTolerance = 1e-4; // relative } - public enum PowellConvergenceType { Success, MaxIterationsExceeded, MaxFunctionCallsExceeded }; + public enum PowellConvergenceType + { + Success, + MaxIterationsExceeded, + MaxFunctionCallsExceeded + } /// /// Result of Powell Minimization. @@ -55,9 +59,9 @@ namespace MathNet.Numerics.Optimization public int NumberOfFunctionCalls; public double[] MinimumPoint; public double MinimumFunctionValue; - public PowellConvergenceType ConvergenceType; + public PowellConvergenceType ConvergenceType; } - + /// /// Minimizes f(p) where p is a vector of model parameters using the Powell method. /// @@ -76,17 +80,17 @@ namespace MathNet.Numerics.Optimization public double[] Minimize(Func function, double[] pInitialGuess) { BrentMinimizer brentMinimizer = new BrentMinimizer(); - int n = pInitialGuess.Length; // number of dimensions + int n = pInitialGuess.Length; // number of dimensions // used in closure: double[] point = new double[n]; double[] startingPoint = new double[n]; double[] direction = new double[n]; double lineMiniumum = 0; int functionCalls = 0; - Func functionAlongLine = (p) => + Func functionAlongLine = p => { for (int i = 0; i < point.Length; ++i) - point[i] = startingPoint[i] + direction[i] * p; + point[i] = startingPoint[i] + direction[i]*p; lineMiniumum = function(point); functionCalls++; return lineMiniumum; @@ -95,8 +99,8 @@ namespace MathNet.Numerics.Optimization double fval; int iterations = 0; - int maxIterations = (Options.MaximumIterations == null) ? n * 1000 : (int)Options.MaximumIterations; - int maxFunctionCalls = (Options.MaximumFunctionCalls == null) ? n * 1000 : (int)Options.MaximumFunctionCalls; + 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[][] directionSet = new double[n][]; @@ -110,7 +114,7 @@ namespace MathNet.Numerics.Optimization double[] x2 = new double[n]; double[] direction1 = new double[n]; - brentMinimizer.Options.Tolerance = Options.PointTolerance * 100; + brentMinimizer.Options.Tolerance = Options.PointTolerance*100; fval = function(x); @@ -142,26 +146,26 @@ namespace MathNet.Numerics.Optimization } } iterations++; - if (2.0 * (fx - fval) <= Options.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 (iterations >= maxIterations) break; - // Construct the extrapolated point + // Construct the extrapolated point for (int i = 0; i < n; ++i) { direction1[i] = x[i] - x1[i]; - x2[i] = 2.0 * x[i] - x1[i]; + x2[i] = 2.0*x[i] - x1[i]; x1[i] = x[i]; } fx2 = function(x2); if (fx > fx2) { - double t = 2.0 * (fx + fx2 - 2.0 * fval); + double t = 2.0*(fx + fx2 - 2.0*fval); double temp = (fx - fval - delta); - t *= temp * temp; + t *= temp*temp; temp = fx - fx2; - t -= delta * temp * temp; + t -= delta*temp*temp; if (t < 0.0) { // Do a linesearch along direction @@ -180,12 +184,12 @@ namespace MathNet.Numerics.Optimization } } var convergenceType = PowellConvergenceType.Success; - if (functionCalls >= maxFunctionCalls) + if (functionCalls >= maxFunctionCalls) convergenceType = PowellConvergenceType.MaxFunctionCallsExceeded; - else if (iterations > maxIterations) + else if (iterations > maxIterations) convergenceType = PowellConvergenceType.MaxFunctionCallsExceeded; - Result = new PowellResult() + Result = new PowellResult { MinimumPoint = (double[])x.Clone(), MinimumFunctionValue = fx, @@ -193,7 +197,7 @@ namespace MathNet.Numerics.Optimization NumberOfIterations = iterations, NumberOfFunctionCalls = functionCalls }; - + return Result.MinimumPoint; } } diff --git a/src/Numerics/Providers/Optimization/IOptimizationProvider.cs b/src/Numerics/Providers/Optimization/IOptimizationProvider.cs index fc137554..4ba98a08 100644 --- a/src/Numerics/Providers/Optimization/IOptimizationProvider.cs +++ b/src/Numerics/Providers/Optimization/IOptimizationProvider.cs @@ -3,9 +3,9 @@ // 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 @@ -14,10 +14,10 @@ // 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 @@ -48,7 +48,7 @@ namespace MathNet.Numerics.Providers.Optimization /// THe model parameters. The function must not change these. /// THe Jacobian matrix in column-major format. public delegate void Jacobian(double[] p, double[] jacobian); - + /// /// Interface to linear algebra algorithms that work off 1-D arrays. /// diff --git a/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs b/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs index ed9a9d84..82510e68 100644 --- a/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs +++ b/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs @@ -172,7 +172,7 @@ namespace MathNet.Numerics.Providers.Optimization.Mkl } SafeNativeMethods.FreeBuffers(); - + NonLinearLeastSquaresConvergenceType convergenceType = NonLinearLeastSquaresConvergenceType.Error; switch (rciRequest) { diff --git a/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs b/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs index 6e4b0735..ad52a84d 100644 --- a/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs +++ b/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs @@ -3,9 +3,9 @@ // 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 @@ -14,10 +14,10 @@ // 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 @@ -28,12 +28,13 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; +using System.Diagnostics; +using MathNet.Numerics.Optimization; +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.OptimizationTests { - using System; - using MathNet.Numerics.Optimization; - using NUnit.Framework; - [TestFixture] public class FunctionMinimizationTests { @@ -44,11 +45,11 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests 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)); + Func function = (x, p) => p[0]*(1 - Math.Exp(-p[1]*x)); var minimizer = new PowellMinimizer(); - var watch = new System.Diagnostics.Stopwatch(); watch.Start(); + var watch = Stopwatch.StartNew(); double[] popt = null; minimizer.Options.PointTolerance = 1e-8; @@ -56,13 +57,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { 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)); + 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-4); Assert.AreEqual(expected[1], popt[1], 1e-4); diff --git a/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs b/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs index 1d458066..1c773b72 100644 --- a/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs +++ b/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs @@ -3,9 +3,9 @@ // 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 @@ -14,10 +14,10 @@ // 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 @@ -30,12 +30,12 @@ #if NATIVEMKL +using System; +using MathNet.Numerics.Optimization; +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.OptimizationTests { - using System; - using MathNet.Numerics.Optimization; - using NUnit.Framework; - [TestFixture] public class NonLinearLeastSquaresTest { @@ -48,13 +48,13 @@ 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 }; - + // 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), + Func jacobian = (x, p) => new double[] { + 1 - Math.Exp(-p[1] * x), p[0] * x * Math.Exp(-p[1] * x) }; var popt = minimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }, jacobian); // 100, 0.75