From 67db4f0501878bc9046a0a5e41559e9f015c98de Mon Sep 17 00:00:00 2001 From: diluculo Date: Wed, 2 Jan 2019 17:41:22 +0900 Subject: [PATCH] Added a converter from IObjectiveModel to IObjectiveFunction. --- .../NonLinearCurveFittingTests.cs | 64 +++++++++++++++++-- src/Numerics/Optimization/IObjectiveModel.cs | 2 + src/Numerics/Optimization/ObjectiveModel.cs | 24 +++++++ .../ObjectiveModels/FittingObjectiveModel.cs | 23 ++++++- 4 files changed, 107 insertions(+), 6 deletions(-) diff --git a/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs index f0787922..c353b83b 100644 --- a/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs +++ b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs @@ -386,6 +386,19 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests } } + [Test] + public void Bfgs_FindMinimum_BoxBod_Unconstrained() + { + var obj = ObjectiveModel.FittingFunction(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); + var solver = new BfgsMinimizer(1e-10, 1e-10, 1e-10, 100); + var result = solver.FindMinimum(obj, BoxBodStart2); + + for (int i = 0; i < result.MinimizingPoint.Count; i++) + { + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6); + } + } + // model : Thurber (https://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml) // f(x; b1 ... b7) = (b1 + b2*x + b3*x^2 + b4*x^3) / (1 + b5*x + b6*x^2 + b7*x^3) // derivatives: @@ -456,7 +469,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests private Vector ThurberPstd = new DenseVector(new double[] { 4.6647963344E+00, 3.9571156086E+01, 2.8698696102E+01, 5.5675370270E+00, 3.1333340687E-02, 1.4984928198E-02, 6.5842344623E-03 }); - private Vector ThurberInitialGuess = new DenseVector(new double[] { 1000.0, 1000.0, 400.0, 40.0, 0.7, 0.3, 0.03 }); + private Vector ThurberStart = new DenseVector(new double[] { 1000.0, 1000.0, 400.0, 40.0, 0.7, 0.3, 0.03 }); + private Vector ThurberLowerBound = new DenseVector(new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }); + private Vector ThurberUpperBound = new DenseVector(new double[] { 1E6, 1E6, 1E6, 1E6, 1E6, 1E6, 1E6 }); private Vector ThurberScales = new DenseVector(new double[7] { 1000, 1000, 400, 40, 0.7, 0.3, 0.03 }); [Test] @@ -464,7 +479,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberPrime, ThurberX, ThurberY); var solver = new LevenbergMarquardtMinimizer(); - var result = solver.FindMinimum(obj, ThurberInitialGuess); + var result = solver.FindMinimum(obj, ThurberStart); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -478,7 +493,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); var solver = new LevenbergMarquardtMinimizer(); - var result = solver.FindMinimum(obj, ThurberInitialGuess); + var result = solver.FindMinimum(obj, ThurberStart); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -494,7 +509,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests scales: ThurberScales, accuracyOrder: 6); var solver = new TrustRegionDogLegMinimizer(); - var result = solver.FindMinimum(obj, ThurberInitialGuess); + var result = solver.FindMinimum(obj, ThurberStart); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -510,7 +525,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests scales: ThurberScales, accuracyOrder: 6); var solver = new TrustRegionNewtonCGMinimizer(); - var result = solver.FindMinimum(obj, ThurberInitialGuess); + var result = solver.FindMinimum(obj, ThurberStart); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -518,5 +533,44 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 3); } } + + [Test] + public void Bfgs_FindMinimum_Thurber_Unconstrained() + { + var obj = ObjectiveModel.FittingFunction(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); + var solver = new BfgsMinimizer(1e-10, 1e-10, 1e-10, 1000); + var result = solver.FindMinimum(obj, ThurberStart); + + for (int i = 0; i < result.MinimizingPoint.Count; i++) + { + AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 6); + } + } + + [Test] + public void BfgsB_FindMinimum_Thurber() + { + var obj = ObjectiveModel.FittingFunction(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); + var solver = new BfgsBMinimizer(1e-10, 1e-10, 1e-10, 1000); + var result = solver.FindMinimum(obj, ThurberLowerBound, ThurberUpperBound, ThurberStart); + + for (int i = 0; i < result.MinimizingPoint.Count; i++) + { + AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 6); + } + } + + [Test] + public void LBfgs_FindMinimum_Thurber() + { + var obj = ObjectiveModel.FittingFunction(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); + var solver = new LimitedMemoryBfgsMinimizer(1e-10, 1e-10, 1e-10, 1000); + var result = solver.FindMinimum(obj, ThurberStart); + + for (int i = 0; i < result.MinimizingPoint.Count; i++) + { + AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 6); + } + } } } diff --git a/src/Numerics/Optimization/IObjectiveModel.cs b/src/Numerics/Optimization/IObjectiveModel.cs index 72c15391..2409dffe 100644 --- a/src/Numerics/Optimization/IObjectiveModel.cs +++ b/src/Numerics/Optimization/IObjectiveModel.cs @@ -66,5 +66,7 @@ namespace MathNet.Numerics.Optimization /// Create a new independent copy of this objective function, evaluated at the same point. IObjectiveModel Fork(); + + IObjectiveFunction ToObjectiveFunction(); } } diff --git a/src/Numerics/Optimization/ObjectiveModel.cs b/src/Numerics/Optimization/ObjectiveModel.cs index 631defae..3d786c18 100644 --- a/src/Numerics/Optimization/ObjectiveModel.cs +++ b/src/Numerics/Optimization/ObjectiveModel.cs @@ -36,5 +36,29 @@ namespace MathNet.Numerics.Optimization objective.SetParameters(lowerBound, upperBound, scales, isFixed); return objective; } + + /// + /// Fitting function with a user supplied jacobian for nonlinear least squares regression by the line search algorithm. + /// + public static IObjectiveFunction FittingFunction(Func, double, double> function, Func, double, Vector> derivatives, + Vector observedX, Vector observedY, Vector weight = null) + { + var objective = new FittingObjectiveModel(function, derivatives); + objective.SetObserved(observedX, observedY, weight); + return objective.ToObjectiveFunction(); + } + + /// + /// Fitting function for nonlinear least squares regression by the line search algorithm. + /// The numerical jacobian with accuracy order is used. + /// + public static IObjectiveFunction FittingFunction(Func, double, double> function, + Vector observedX, Vector observedY, Vector weight = null, + int accuracyOrder = 2) + { + var objective = new FittingObjectiveModel(function, null, accuracyOrder: accuracyOrder); + objective.SetObserved(observedX, observedY, weight); + return objective.ToObjectiveFunction(); + } } } diff --git a/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs b/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs index 6b9fba7c..c5558f99 100644 --- a/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs +++ b/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs @@ -1,4 +1,5 @@ using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.Optimization.ObjectiveFunctions; using System; using System.Collections.Generic; using System.Linq; @@ -172,7 +173,27 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels public IObjectiveModel CreateNew() { - return new FittingObjectiveModel(userFunction, userDerivatives); + return new FittingObjectiveModel(userFunction, userDerivatives, AccuracyOrder); + } + + public IObjectiveFunction ToObjectiveFunction() + { + Tuple, Matrix> function(Vector point) + { + EvaluateFunction(point); + EvaluateJacobian(point); + + return new Tuple, Matrix>(Residue, -Gradient, Hessian); + } + + LowerBound = null; + UpperBound = null; + Scales = null; + IsFixed = null; + IsBounded = false; + + var objective = new GradientHessianObjectiveFunction(function); + return objective; } ///