tsai
/
mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Updated the box-constrained mapping rule in the FittingObjectiveModel.

pull/614/head
diluculo 8 years ago
parent
aca34e13b9
commit
220ab5aadb
2 changed files with 155 additions and 75 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 178
      src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs
  2. 52
      src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs

178
src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs
View File

@ -42,11 +42,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
{
var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY);
var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 3);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 3);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 3);
}
}
[Test]
@ -54,11 +55,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
{
var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockX, RosenbrockY, accuracyOrder:2);
var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 3);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 3);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 3);
}
}
[Test]
@ -67,11 +69,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY,
lowerBound : RosebbrockLowerBound, upperBound : RosenbrockUpperBound);
var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 3);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 3);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 3);
}
}
[Test]
@ -79,13 +82,14 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
{
var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockX, RosenbrockY,
lowerBound: RosebbrockLowerBound, upperBound: RosenbrockUpperBound,
accuracyOrder: 2);
accuracyOrder: 6);
var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 3);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 3);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 3);
}
}
[Test]
@ -93,11 +97,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
{
var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY);
var solver = new TrustRegionDogLegMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 1);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 1);
}
}
[Test]
@ -105,11 +110,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
{
var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY);
var solver = new TrustRegionDogLegMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 1);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 1);
}
}
// model: Rat43 (https://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml)
@ -147,7 +153,6 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
{
var obj = ObjectiveModel.FittingModel(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, Rat43Start1);
for (int i = 0; i < result.BestFitParameters.Count; i++)
@ -162,7 +167,6 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
{
var obj = ObjectiveModel.FittingModel(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6);
var solver = new TrustRegionDogLegMinimizer();
var result = solver.FindMinimum(obj, Rat43Start2);
for (int i = 0; i < result.BestFitParameters.Count; i++)
@ -199,23 +203,22 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
private Vector<double> BoxBodStart1 = new DenseVector(new double[] { 1.0, 1.0 });
private Vector<double> BoxBodStart2 = new DenseVector(new double[] { 100.0, 0.75 });
private Vector<double> BoxBodLowerBound = new DenseVector(new double[] { 0, 0 });
private Vector<double> BoxBodUpperBound = new DenseVector(new double[] { 500.0, 10 });
private Vector<double> BoxBodScales = new DenseVector(new double[] { 100.0, 1 });
private Vector<double> BoxBodLowerBound = new DenseVector(new double[] { -1000, -100 });
private Vector<double> BoxBodUpperBound = new DenseVector(new double[] { 1000.0, 100 });
private Vector<double> BoxBodScales = new DenseVector(new double[] { 100.0, 0.1 });
[Test]
public void LMDER_FindMinimum_BoxBod_Unconstrained()
{
var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 6);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
}
[Test]
@ -223,30 +226,96 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
{
var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder:6);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 6);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
}
[Test]
public void LMDER_FindMinimum_BoxBod_BoxConstrained()
{
// lower < parameters < upper
// Note that in this case, scales have no effect.
var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY,
lowerBound: BoxBodLowerBound, upperBound: BoxBodUpperBound);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 6);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 6);
// lower < parameters, no scales
obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY,
lowerBound: BoxBodLowerBound);
solver = new LevenbergMarquardtMinimizer();
result = solver.FindMinimum(obj, BoxBodStart1);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
// lower < parameters, scales
obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY,
lowerBound: BoxBodLowerBound, scales: BoxBodScales);
solver = new LevenbergMarquardtMinimizer();
result = solver.FindMinimum(obj, BoxBodStart1);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
// parameters < upper, no scales
obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY,
upperBound: BoxBodUpperBound);
solver = new LevenbergMarquardtMinimizer();
result = solver.FindMinimum(obj, BoxBodStart1);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
// parameters < upper, scales
obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY,
upperBound: BoxBodUpperBound, scales: BoxBodScales);
solver = new LevenbergMarquardtMinimizer();
result = solver.FindMinimum(obj, BoxBodStart1);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
// only scales
obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY,
scales: BoxBodScales);
solver = new LevenbergMarquardtMinimizer();
result = solver.FindMinimum(obj, BoxBodStart1);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
}
[Test]
@ -256,14 +325,13 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
lowerBound: BoxBodLowerBound, upperBound: BoxBodUpperBound,
accuracyOrder: 6);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 6);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
}
[Test]
@ -271,14 +339,13 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
{
var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6);
var solver = new TrustRegionDogLegMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 3);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 3);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 3);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 3);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 3);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 3);
}
}
// model : Thurber (https://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml)
@ -359,7 +426,6 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
{
var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberPrime, ThurberX, ThurberY);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, ThurberInitialGuess);
for (int i = 0; i < result.BestFitParameters.Count; i++)
@ -374,7 +440,6 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
{
var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, ThurberInitialGuess);
for (int i = 0; i < result.BestFitParameters.Count; i++)
@ -391,7 +456,6 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
scales: ThurberScales,
accuracyOrder: 6);
var solver = new TrustRegionDogLegMinimizer();
var result = solver.FindMinimum(obj, ThurberInitialGuess);
for (int i = 0; i < result.BestFitParameters.Count; i++)

52
src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs
View File

@ -329,30 +329,34 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels
ValidateParameters(parameters);
// To handle the box constrained minimization as the unconstrained minimization,
// parameters are mapping by the following rule.
// the parameters are mapping by the following rules,
// which are modified the rules shown in the ref[1] in order to introduce scales.
//
// 1. lower < P < upper
// 1. lower < Pext < upper
// Pint = asin(2 * (Pext - lower) / (upper - lower) - 1)
// Pext = lower + (sin(Pint) + 1) * (upper - lower) / 2
// dPext/dPint = (upper - lower) / 2 * cos(Pint)
// 2. lower < P
// Pint = sqrt((Pext - lower + 1)^2 - 1)
// Pext = lower - 1 + sqrt(Pint^2 + 1)
// dPext/dPint = Pint / sqrt(Pint^2 + 1)
// 3. P < upper
// Pint = sqrt((upper - Pext + 1)^2 - 1)
// Pext = upper + 1 - sqrt(Pint^2 + 1)
// dPext/dPint = - Pint / sqrt(Pint^2 + 1)
//
// 2. lower < Pext
// Pint = sqrt((Pext/scale - lower/scale + 1)^2 - 1)
// Pext = lower + scale * (sqrt(Pint^2 + 1) - 1)
// dPext/dPint = scale * Pint / sqrt(Pint^2 + 1)
//
// 3. Pext < upper
// Pint = sqrt((upper / scale - Pext / scale + 1)^2 - 1)
// Pext = upper + scale - scale * sqrt(Pint^2 + 1)
// dPext/dPint = - scale * Pint / sqrt(Pint^2 + 1)
//
// 4. no bounds, but scales
// Pint = Pext / scale
// Pext = Pint * scale
// dPext/dPint = scale
//
// see ProjectParametersToInternal(Pext), ProjectParametersToExternal(Pint), ScaleFactorsOfJacobian(Pint) methods.
// The rules are applied in ProjectParametersToInternal, ProjectParametersToExternal, and ScaleFactorsOfJacobian methods.
//
// References:
// [1] https://lmfit.github.io/lmfit-py/bounds.html
//
// [2] MINUIT User's Guide, https://root.cern.ch/download/minuit.pdf
//
// Except when it is initial guess, the parameters argument is always internal parameter.
// So, first map the parameters argument to the external parameters in order to calculate function values.
@ -617,7 +621,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels
{
for (int i = 0; i < Pext.Count; i++)
{
Pint[i] = Math.Sqrt(Math.Pow(Pext[i] - LowerBound[i] + 1.0, 2) - 1.0);
Pint[i] = (Scales == null)
? Math.Sqrt(Math.Pow(Pext[i] - LowerBound[i] + 1.0, 2) - 1.0)
: Math.Sqrt(Math.Pow((Pext[i] - LowerBound[i]) / Scales[i] + 1.0, 2) - 1.0);
}
return Pint;
@ -626,7 +632,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels
{
for (int i = 0; i < Pext.Count; i++)
{
Pint[i] = Math.Sqrt(Math.Pow(UpperBound[i] - Pext[i] + 1.0, 2) - 1.0);
Pint[i] = (Scales == null)
? Math.Sqrt(Math.Pow(UpperBound[i] - Pext[i] + 1.0, 2) - 1.0)
: Math.Sqrt(Math.Pow((UpperBound[i] - Pext[i]) / Scales[i] + 1.0, 2) - 1.0);
}
return Pint;
@ -661,7 +669,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels
{
for (int i = 0; i < Pint.Count; i++)
{
Pext[i] = LowerBound[i] + Math.Sqrt(Pint[i] * Pint[i] + 1.0) - 1.0;
Pext[i] = (Scales == null)
? LowerBound[i] + Math.Sqrt(Pint[i] * Pint[i] + 1.0) - 1.0
: LowerBound[i] + Scales[i] * (Math.Sqrt(Pint[i] * Pint[i] + 1.0) - 1.0);
}
return Pext;
@ -670,7 +680,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels
{
for (int i = 0; i < Pint.Count; i++)
{
Pext[i] = UpperBound[i] - Math.Sqrt(Pint[i] * Pint[i] + 1.0) + 1.0;
Pext[i] = (Scales == null)
? UpperBound[i] - Math.Sqrt(Pint[i] * Pint[i] + 1.0) + 1.0
: UpperBound[i] - Scales[i] * (Math.Sqrt(Pint[i] * Pint[i] + 1.0) - 1.0);
}
return Pext;
@ -704,7 +716,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels
{
for (int i = 0; i < Pint.Count; i++)
{
scale[i] = Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0);
scale[i] = (Scales == null)
? Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0)
: Scales[i] * Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0);
}
return scale;
}
@ -712,7 +726,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels
{
for (int i = 0; i < Pint.Count; i++)
{
scale[i] = -Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0);
scale[i] = (Scales == null)
? -Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0)
: -Scales[i] * Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0);
}
return scale;
}

Write Preview
Loading…
Cancel
Save