mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-01 18:26:24 +08:00
Code Issues Packages Projects Releases Wiki Activity

NonLinearOptimization : clean 'mode' handling from the old minpack code :

Browse Source 
* this is actually a boolean, not an int
* use a better name
* can be set at initialization time instead of bloating all methods signatures
This commit is contained in:
Thomas Capricelli 2010-02-21 12:41:37 +01:00
parent 608959aa6f
commit a7d085eb4e
3 changed files with 68 additions and 141 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

108
unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
View File

@ -57,7 +57,7 @@ class HybridNonLinearSolver
{
public:
HybridNonLinearSolver(FunctorType &_functor)
: functor(_functor) { nfev=njev=iter = 0; fnorm= 0.; }
: functor(_functor) { nfev=njev=iter = 0; fnorm= 0.; useExternalScaling=false;}
struct Parameters {
Parameters()
@ -84,36 +84,18 @@ public:
const Scalar tol = ei_sqrt(NumTraits<Scalar>::epsilon())
);
HybridNonLinearSolverSpace::Status solveInit(
FVectorType &x,
const int mode=1
);
HybridNonLinearSolverSpace::Status solveOneStep(
FVectorType &x,
const int mode=1
);
HybridNonLinearSolverSpace::Status solve(
FVectorType &x,
const int mode=1
);
HybridNonLinearSolverSpace::Status solveInit(FVectorType &x);
HybridNonLinearSolverSpace::Status solveOneStep(FVectorType &x);
HybridNonLinearSolverSpace::Status solve(FVectorType &x);
HybridNonLinearSolverSpace::Status hybrd1(
FVectorType &x,
const Scalar tol = ei_sqrt(NumTraits<Scalar>::epsilon())
);
HybridNonLinearSolverSpace::Status solveNumericalDiffInit(
FVectorType &x,
const int mode=1
);
HybridNonLinearSolverSpace::Status solveNumericalDiffOneStep(
FVectorType &x,
const int mode=1
);
HybridNonLinearSolverSpace::Status solveNumericalDiff(
FVectorType &x,
const int mode=1
);
HybridNonLinearSolverSpace::Status solveNumericalDiffInit(FVectorType &x);
HybridNonLinearSolverSpace::Status solveNumericalDiffOneStep(FVectorType &x);
HybridNonLinearSolverSpace::Status solveNumericalDiff(FVectorType &x);
void resetParameters(void) { parameters = Parameters(); }
Parameters parameters;
@ -124,6 +106,7 @@ public:
int njev;
int iter;
Scalar fnorm;
bool useExternalScaling;
private:
FunctorType &functor;
int n;
@ -160,18 +143,13 @@ HybridNonLinearSolver<FunctorType,Scalar>::hybrj1(
parameters.maxfev = 100*(n+1);
parameters.xtol = tol;
diag.setConstant(n, 1.);
return solve(
x,
2
);
useExternalScaling = true;
return solve(x);
}
template<typename FunctorType, typename Scalar>
HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
FVectorType &x,
const int mode
)
HybridNonLinearSolver<FunctorType,Scalar>::solveInit(FVectorType &x)
{
n = x.size();
@ -179,9 +157,9 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
fvec.resize(n);
qtf.resize(n);
fjac.resize(n, n);
if (mode != 2)
if (!useExternalScaling)
diag.resize(n);
assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'");
assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
/* Function Body */
nfev = 0;
@ -190,7 +168,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
/* check the input parameters for errors. */
if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0. )
return HybridNonLinearSolverSpace::ImproperInputParameters;
if (mode == 2)
if (useExternalScaling)
for (int j = 0; j < n; ++j)
if (diag[j] <= 0.)
return HybridNonLinearSolverSpace::ImproperInputParameters;
@ -214,10 +192,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
template<typename FunctorType, typename Scalar>
HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
FVectorType &x,
const int mode
)
HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(FVectorType &x)
{
int j;
std::vector<PlanarRotation<Scalar> > v_givens(n), w_givens(n);
@ -231,10 +206,10 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
wa2 = fjac.colwise().blueNorm();
/* on the first iteration and if mode is 1, scale according */
/* on the first iteration and if external scaling is not used, scale according */
/* to the norms of the columns of the initial jacobian. */
if (iter == 1) {
if (mode != 2)
if (!useExternalScaling)
for (j = 0; j < n; ++j)
diag[j] = (wa2[j]==0.) ? 1. : wa2[j];
@ -260,7 +235,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
qtf = fjac.transpose() * fvec;
/* rescale if necessary. */
if (mode != 2)
if (!useExternalScaling)
diag = diag.cwiseMax(wa2);
while (true) {
@ -372,14 +347,11 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
template<typename FunctorType, typename Scalar>
HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solve(
FVectorType &x,
const int mode
)
HybridNonLinearSolver<FunctorType,Scalar>::solve(FVectorType &x)
{
HybridNonLinearSolverSpace::Status status = solveInit(x, mode);
HybridNonLinearSolverSpace::Status status = solveInit(x);
while (status==HybridNonLinearSolverSpace::Running)
status = solveOneStep(x, mode);
status = solveOneStep(x);
return status;
}
@ -403,18 +375,13 @@ HybridNonLinearSolver<FunctorType,Scalar>::hybrd1(
parameters.xtol = tol;
diag.setConstant(n, 1.);
return solveNumericalDiff(
x,
2
);
useExternalScaling = true;
return solveNumericalDiff(x);
}
template<typename FunctorType, typename Scalar>
HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
FVectorType &x,
const int mode
)
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(FVectorType &x)
{
n = x.size();
@ -425,10 +392,9 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
qtf.resize(n);
fjac.resize(n, n);
fvec.resize(n);
if (mode != 2)
if (!useExternalScaling)
diag.resize(n);
assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'");
assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
/* Function Body */
nfev = 0;
@ -437,7 +403,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
/* check the input parameters for errors. */
if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.nb_of_subdiagonals< 0 || parameters.nb_of_superdiagonals< 0 || parameters.factor <= 0. )
return HybridNonLinearSolverSpace::ImproperInputParameters;
if (mode == 2)
if (useExternalScaling)
for (int j = 0; j < n; ++j)
if (diag[j] <= 0.)
return HybridNonLinearSolverSpace::ImproperInputParameters;
@ -461,10 +427,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
template<typename FunctorType, typename Scalar>
HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
FVectorType &x,
const int mode
)
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(FVectorType &x)
{
int j;
std::vector<PlanarRotation<Scalar> > v_givens(n), w_givens(n);
@ -480,10 +443,10 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
wa2 = fjac.colwise().blueNorm();
/* on the first iteration and if mode is 1, scale according */
/* on the first iteration and if external scaling is not used, scale according */
/* to the norms of the columns of the initial jacobian. */
if (iter == 1) {
if (mode != 2)
if (!useExternalScaling)
for (j = 0; j < n; ++j)
diag[j] = (wa2[j]==0.) ? 1. : wa2[j];
@ -509,7 +472,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
qtf = fjac.transpose() * fvec;
/* rescale if necessary. */
if (mode != 2)
if (!useExternalScaling)
diag = diag.cwiseMax(wa2);
while (true) {
@ -621,14 +584,11 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
template<typename FunctorType, typename Scalar>
HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiff(
FVectorType &x,
const int mode
)
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiff(FVectorType &x)
{
HybridNonLinearSolverSpace::Status status = solveNumericalDiffInit(x, mode);
HybridNonLinearSolverSpace::Status status = solveNumericalDiffInit(x);
while (status==HybridNonLinearSolverSpace::Running)
status = solveNumericalDiffOneStep(x, mode);
status = solveNumericalDiffOneStep(x);
return status;
}

95
unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
View File

@ -61,7 +61,7 @@ class LevenbergMarquardt
{
public:
LevenbergMarquardt(FunctorType &_functor)
: functor(_functor) { nfev = njev = iter = 0; fnorm=gnorm = 0.; }
: functor(_functor) { nfev = njev = iter = 0; fnorm = gnorm = 0.; useExternalScaling=false; }
struct Parameters {
Parameters()
@ -87,18 +87,9 @@ public:
const Scalar tol = ei_sqrt(NumTraits<Scalar>::epsilon())
);
LevenbergMarquardtSpace::Status minimize(
FVectorType &x,
const int mode=1
);
LevenbergMarquardtSpace::Status minimizeInit(
FVectorType &x,
const int mode=1
);
LevenbergMarquardtSpace::Status minimizeOneStep(
FVectorType &x,
const int mode=1
);
LevenbergMarquardtSpace::Status minimize(FVectorType &x);
LevenbergMarquardtSpace::Status minimizeInit(FVectorType &x);
LevenbergMarquardtSpace::Status minimizeOneStep(FVectorType &x);
static LevenbergMarquardtSpace::Status lmdif1(
FunctorType &functor,
@ -112,18 +103,9 @@ public:
const Scalar tol = ei_sqrt(NumTraits<Scalar>::epsilon())
);
LevenbergMarquardtSpace::Status minimizeOptimumStorage(
FVectorType &x,
const int mode=1
);
LevenbergMarquardtSpace::Status minimizeOptimumStorageInit(
FVectorType &x,
const int mode=1
);
LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep(
FVectorType &x,
const int mode=1
);
LevenbergMarquardtSpace::Status minimizeOptimumStorage(FVectorType &x);
LevenbergMarquardtSpace::Status minimizeOptimumStorageInit(FVectorType &x);
LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep(FVectorType &x);
void resetParameters(void) { parameters = Parameters(); }
@ -135,6 +117,7 @@ public:
int njev;
int iter;
Scalar fnorm, gnorm;
bool useExternalScaling;
Scalar lm_param(void) { return par; }
private:
@ -175,24 +158,18 @@ LevenbergMarquardt<FunctorType,Scalar>::lmder1(
template<typename FunctorType, typename Scalar>
LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimize(
FVectorType &x,
const int mode
)
LevenbergMarquardt<FunctorType,Scalar>::minimize(FVectorType &x)
{
LevenbergMarquardtSpace::Status status = minimizeInit(x, mode);
LevenbergMarquardtSpace::Status status = minimizeInit(x);
do {
status = minimizeOneStep(x, mode);
status = minimizeOneStep(x);
} while (status==LevenbergMarquardtSpace::Running);
return status;
}
template<typename FunctorType, typename Scalar>
LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
FVectorType &x,
const int mode
)
LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(FVectorType &x)
{
n = x.size();
m = functor.values();
@ -201,9 +178,9 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
wa4.resize(m);
fvec.resize(m);
fjac.resize(m, n);
if (mode != 2)
if (!useExternalScaling)
diag.resize(n);
assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'");
assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
qtf.resize(n);
/* Function Body */
@ -214,7 +191,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
return LevenbergMarquardtSpace::ImproperInputParameters;
if (mode == 2)
if (useExternalScaling)
for (int j = 0; j < n; ++j)
if (diag[j] <= 0.)
return LevenbergMarquardtSpace::ImproperInputParameters;
@ -235,10 +212,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
template<typename FunctorType, typename Scalar>
LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
FVectorType &x,
const int mode
)
LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(FVectorType &x)
{
int j;
@ -257,10 +231,10 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
fjac = qrfac.matrixQR();
permutation = qrfac.colsPermutation();
/* on the first iteration and if mode is 1, scale according */
/* on the first iteration and if external scaling is not used, scale according */
/* to the norms of the columns of the initial jacobian. */
if (iter == 1) {
if (mode != 2)
if (!useExternalScaling)
for (j = 0; j < n; ++j)
diag[j] = (wa2[j]==0.)? 1. : wa2[j];
@ -290,7 +264,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
return LevenbergMarquardtSpace::CosinusTooSmall;
/* rescale if necessary. */
if (mode != 2)
if (!useExternalScaling)
diag = diag.cwiseMax(wa2);
do {
@ -406,10 +380,7 @@ LevenbergMarquardt<FunctorType,Scalar>::lmstr1(
template<typename FunctorType, typename Scalar>
LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
FVectorType &x,
const int mode
)
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(FVectorType &x)
{
n = x.size();
m = functor.values();
@ -423,9 +394,9 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
// The purpose it to only use a nxn matrix, instead of mxn here, so
// that we can handle cases where m>>n :
fjac.resize(n, n);
if (mode != 2)
if (!useExternalScaling)
diag.resize(n);
assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'");
assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
qtf.resize(n);
/* Function Body */
@ -436,7 +407,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
return LevenbergMarquardtSpace::ImproperInputParameters;
if (mode == 2)
if (useExternalScaling)
for (int j = 0; j < n; ++j)
if (diag[j] <= 0.)
return LevenbergMarquardtSpace::ImproperInputParameters;
@ -458,10 +429,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
template<typename FunctorType, typename Scalar>
LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
FVectorType &x,
const int mode
)
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(FVectorType &x)
{
int i, j;
bool sing;
@ -514,10 +482,10 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
}
}
/* on the first iteration and if mode is 1, scale according */
/* on the first iteration and if external scaling is not used, scale according */
/* to the norms of the columns of the initial jacobian. */
if (iter == 1) {
if (mode != 2)
if (!useExternalScaling)
for (j = 0; j < n; ++j)
diag[j] = (wa2[j]==0.)? 1. : wa2[j];
@ -541,7 +509,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
return LevenbergMarquardtSpace::CosinusTooSmall;
/* rescale if necessary. */
if (mode != 2)
if (!useExternalScaling)
diag = diag.cwiseMax(wa2);
do {
@ -635,14 +603,11 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
template<typename FunctorType, typename Scalar>
LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(
FVectorType &x,
const int mode
)
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(FVectorType &x)
{
LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x, mode);
LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x);
do {
status = minimizeOptimumStorageOneStep(x, mode);
status = minimizeOptimumStorageOneStep(x);
} while (status==LevenbergMarquardtSpace::Running);
return status;
}

6
unsupported/test/NonLinearOptimization.cpp
View File

@ -317,7 +317,8 @@ void testHybrj()
hybrj_functor functor;
HybridNonLinearSolver<hybrj_functor> solver(functor);
solver.diag.setConstant(n, 1.);
info = solver.solve(x, 2);
solver.useExternalScaling = true;
info = solver.solve(x);
// check return value
VERIFY( 1 == info);
@ -401,7 +402,8 @@ void testHybrd()
solver.parameters.nb_of_subdiagonals = 1;
solver.parameters.nb_of_superdiagonals = 1;
solver.diag.setConstant(n, 1.);
info = solver.solveNumericalDiff(x, 2);
solver.useExternalScaling = true;
info = solver.solveNumericalDiff(x);
// check return value
VERIFY( 1 == info);