diff --git a/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs b/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs index 5611477c..bec4f902 100644 --- a/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs +++ b/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs @@ -6,66 +6,77 @@ using MathNet.Numerics.Providers.Optimization; using MathNet.Numerics.Providers.Optimization.Mkl; 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. + /// Options for Non-Linear Least Squares Minimization. /// - public class NonLinearLeastSquaresMinimizer + public class NonLinearLeastSquaresOptions { - public class Options - { - public int MaximumIterations = 1000; + public int MaximumIterations = 1000; - public int MaximumTrialStepIterations = 100; + public int MaximumTrialStepIterations = 100; - public ConvergenceType ConvergenceType; + public NonLinearLeastSquaresConvergenceType ConvergenceType; - /// - /// Convergence if Δ < Criterion0, Δ is trust region size. - /// - public double Criterion0 = 1e-7; + /// + /// Convergence if Δ < Criterion0, Δ is trust region size. + /// + public double Criterion0 = 1e-7; - /// - /// Convergence if |r|2 < Criterion1, r is residuals vector. - /// - public double Criterion1 = 1e-7; + /// + /// Convergence if |r|2 < Criterion1, r is residuals vector. + /// + public double Criterion1 = 1e-7; - /// - /// Jacobian considered singular if |J(:,j)|2 < Criterion2 for any j. - /// - public double Criterion2 = 1e-7; + /// + /// Jacobian considered singular if |J(:,j)|2 < Criterion2 for any j. + /// + public double Criterion2 = 1e-7; - /// - /// Jacobian considered singular if |s|2 < Criterion3. s is trial step. - /// - public double Criterion3 = 1e-7; + /// + /// Jacobian considered singular if |s|2 < Criterion3. s is trial step. + /// + public double Criterion3 = 1e-7; - /// - /// |r|2 - |r - Js|2 < Criterion4 - /// - public double Criterion4 = 1e-7; + /// + /// |r|2 - |r - Js|2 < Criterion4 + /// + public double Criterion4 = 1e-7; - public double TrialStepPrecision = 1e-10; + public double TrialStepPrecision = 1e-10; - /// - /// Only if Jacobian calculated by central differences. - /// - public double JacobianPrecision = 1e-8; - } - /// - /// For details of convergence criteria, see Options. + /// Only if Jacobian calculated by central differences. /// - public enum ConvergenceType { NoneMaxIterationExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, Error }; - - public class Result - { - public int NumberOfIterations; + public double JacobianPrecision = 1e-8; + } + + /// + /// For details of convergence criteria, see Options. + /// + public enum NonLinearLeastSquaresConvergenceType { NoneMaxIterationExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, Error }; + + /// + /// Result of Non-Linear Least Squares Minimization. + /// + public class NonLinearLeastSquaresResult + { + public int NumberOfIterations; + + public NonLinearLeastSquaresConvergenceType ConvergenceType; + } + + + /// + /// 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. + /// + public class NonLinearLeastSquaresMinimizer + { + public NonLinearLeastSquaresResult Result { get; private set; } + + public readonly NonLinearLeastSquaresOptions Options = new NonLinearLeastSquaresOptions(); - public ConvergenceType ConvergenceType; - } - /// /// Non-Linear Least-Squares fitting the points (x,y) to a specified function of y : x -> f(x, p), p being a vector of parameters. /// returning its best fitting parameters p. @@ -76,7 +87,7 @@ namespace MathNet.Numerics.Optimization /// Initial guess of parameters, p. /// jac_j(x, p) = df / dp_j /// - public static double[] CurveFit(double[] x, double[] y, Func f, + public double[] CurveFit(double[] x, double[] y, Func f, double[] pStart, Func jacobian = null) { if (x.Length != y.Length) throw new ArgumentException("x and y lengths different"); @@ -101,7 +112,7 @@ namespace MathNet.Numerics.Optimization }; double[] parameters; - Result result = provider.NonLinearLeastSquaresUnboundedMinimize(y.Length, pStart, function, out parameters, jacobianFunction); + Result = provider.NonLinearLeastSquaresUnboundedMinimize(y.Length, pStart, function, out parameters, jacobianFunction); return parameters; } } diff --git a/src/Numerics/Providers/Optimization/IOptimizationProvider.cs b/src/Numerics/Providers/Optimization/IOptimizationProvider.cs index 3e2400ba..fc137554 100644 --- a/src/Numerics/Providers/Optimization/IOptimizationProvider.cs +++ b/src/Numerics/Providers/Optimization/IOptimizationProvider.cs @@ -56,8 +56,8 @@ namespace MathNet.Numerics.Providers.Optimization public interface IOptimizationProvider where T : struct { - NonLinearLeastSquaresMinimizer.Result NonLinearLeastSquaresUnboundedMinimize( + NonLinearLeastSquaresResult NonLinearLeastSquaresUnboundedMinimize( int residualsLength, T[] initialGuess, LeastSquaresForwardModel function, - out T[] parameters, Jacobian jacobianFunction = null, NonLinearLeastSquaresMinimizer.Options options = null); + out T[] parameters, Jacobian jacobianFunction = null, NonLinearLeastSquaresOptions options = null); } } diff --git a/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs b/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs index eeac095e..b5530720 100644 --- a/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs +++ b/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs @@ -42,9 +42,9 @@ namespace MathNet.Numerics.Providers.Optimization.Mkl { const int TR_SUCCESS = 1501; - public NonLinearLeastSquaresMinimizer.Result NonLinearLeastSquaresUnboundedMinimize(int residualsLength, double[] initialGuess, LeastSquaresForwardModel function, out double[] parameters, Jacobian jacobianFunction = null, NonLinearLeastSquaresMinimizer.Options options = null) + public NonLinearLeastSquaresResult NonLinearLeastSquaresUnboundedMinimize(int residualsLength, double[] initialGuess, LeastSquaresForwardModel function, out double[] parameters, Jacobian jacobianFunction = null, NonLinearLeastSquaresOptions options = null) { - if (options == null) options = new NonLinearLeastSquaresMinimizer.Options(); + if (options == null) options = new NonLinearLeastSquaresOptions(); bool analyticJacobian = jacobianFunction != null; double[] residuals = new double[residualsLength]; double[] residualsMinus = new double[residualsLength]; @@ -176,30 +176,30 @@ namespace MathNet.Numerics.Providers.Optimization.Mkl SafeNativeMethods.FreeBuffers(); - NonLinearLeastSquaresMinimizer.ConvergenceType convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Error; + NonLinearLeastSquaresConvergenceType convergenceType = NonLinearLeastSquaresConvergenceType.Error; switch (rciRequest) { case -1: - convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.NoneMaxIterationExceeded; break; + convergenceType = NonLinearLeastSquaresConvergenceType.NoneMaxIterationExceeded; break; case -2: - convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion0; break; + convergenceType = NonLinearLeastSquaresConvergenceType.Criterion0; break; case -3: - convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion1; break; + convergenceType = NonLinearLeastSquaresConvergenceType.Criterion1; break; case -4: - convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion2; break; + convergenceType = NonLinearLeastSquaresConvergenceType.Criterion2; break; case -5: - convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion3; break; + convergenceType = NonLinearLeastSquaresConvergenceType.Criterion3; break; case -6: - convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion4; break; + convergenceType = NonLinearLeastSquaresConvergenceType.Criterion4; break; } // no errors, find reason for stopping; - return new NonLinearLeastSquaresMinimizer.Result() { ConvergenceType = convergenceType, NumberOfIterations = iterations }; + return new NonLinearLeastSquaresResult() { ConvergenceType = convergenceType, NumberOfIterations = iterations }; } - public static NonLinearLeastSquaresMinimizer.Result ErrorResult() + public static NonLinearLeastSquaresResult ErrorResult() { - return new NonLinearLeastSquaresMinimizer.Result() { ConvergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Error }; + return new NonLinearLeastSquaresResult() { ConvergenceType = NonLinearLeastSquaresConvergenceType.Error }; } } } diff --git a/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs b/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs index 17cf12d7..325275ef 100644 --- a/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs +++ b/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs @@ -40,17 +40,20 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void CurveFit() { + var minimizer = new NonLinearLeastSquaresMinimizer(); + minimizer.Options.MaximumIterations = 1010; + minimizer.Options.Criterion4 = 1e-8; // 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 = NonLinearLeastSquaresMinimizer.CurveFit(xin, yin, (x, p) => p[0] * (1 - Math.Exp(-p[1] * x)), new double[] { 1, 1 }); + 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 = NonLinearLeastSquaresMinimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }, jacobian); // 100, 0.75 + popt = minimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }, jacobian); // 100, 0.75 double[] expected = new double[] { 2.1380940889E+02, 5.4723748542E-01 };